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Preface 


This  thesis  evolved  from  a  problem  suggested  by  the  Rome 
Air  Development  Center.  The  original  problem  considered  Max¬ 
imum  Entropy  spectral  analysis  of  Non-White  Atmospheric  Radio 
Noise.  RADC  suggested  an  extension  of  the  analysis  to  include 
the  work  by  Papoulis  in  bandlimited  extrapolation.  When  actual 
noise  tapes  were  not  received  from  RADC,  I  changed  the  emphasis 
to  analysis  of  short  record  data.  Data  that  is  extracted  from 
only  a  few  input  samples.  Specific  problems  were  extracted 
from  the  proceedings  of  a  Spectral  Estimation  Workshop  at  RADC 
in  1978.  This  report  covers  the  areas  of  research  to  analyze 
the  MEM  and  Papoulis  spectral  estimation  techniques. 

During  the  course  of  this  work,  valuable  guidance,  in¬ 
sight,  and  probing  questions  were  provided  by  my  thesis  ad¬ 
visor,  Lt.  Col.  Ronald  J.  Carpinella.  His  help  is  most  grate¬ 
fully  acknowledged.  Special  thanks  are  also  due  to  the  read¬ 
ers  Dr.  Pater  S.  Maybeck  and  Capt.  Stanley  R.  Robinson  of  the 
Electrical  Engineering  Department  of  AFIT. 

This  thesis  could  not  have  been  possible  without  the  help 
and  encouragement  of  my  wife  Shauna.  Her  support  and  unself¬ 
ishness  have  enabled  me  to  achieve  a  major  goal  I  have  worked 
for  for  so  many  years. 
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Abstract 


Spectral  estimation  of  data  from  some  radar  applications 
and  seismilogical  events  is  not  accurate  when  short  records 
are  evaluated  using  traditional  techniques.  A  record  of  data 
is  short  if  the  number  of  samples  from  the  process  is  more 
than  an  order  of  magnitude  smaller  than  the  reciprical  of  the 
lowest  frequency  of  interest.  This  analysis  considers  records 
of  fewer  than  128  samples. 

Techniques  that  produce  improved  frequency  and  ampli¬ 
tude  resolution  over  smoothed  periodograms  and  Fast  Fourier 
Transforms  (FFT)  are  considered.  Specifically,  the  Burg  Max¬ 
imum  Entropy  Method  (MEM)  and  Papoulis  Bandlimited  Extrapo¬ 
lation  are  derived.  These  techniques  are  shown  to  produce 
estimated  that  become  unbiased  and  consistant.  Additionally, 
the  effects  of  windowing,  a  problem  inherent  with  periodo¬ 
grams,  are  not  observed  in  these  techniques.. 

Model  order  determination  for  the  MEM  technique  is  accom¬ 


plished  through  implementation  of  two  techniques  developed  by 


Akaike.  These  techniques  are  derived  through  the  maximum  like- 
lihood  estimation  of  the  estimated  power".  Using  Akaike’ s  or¬ 
der  techniques,  MEM  spectral  estimation  provides  extremely 
good  frequency  resolution  from  very  short  input  records. 

^-^Papoulis  bandlimited  extrapolation  techniques  provided 
accurate  results  when  short  records  are  evaluated.  The 


x 


technique  is  limited,  in  this  research,  to  computational  times 
commensurate  with  the  MEM  technique  and  does  not  produce  resol¬ 
ution  accuracy  that  is  observed  from  MEM  analysis.  However, 
the  estimated  amplitudes  are  much  closer  than  those  generated 
from  the  MEM  technique. 
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Spectral  Analysis  of  Short  Record 
Time  Series  Data 

I .  Introduction 

Considerable  interest  is  generated  by  government  research 
agencies  in  power  spectral  estimation  of  time  series  data. 
Traditional  techniques  (Fast  Fourier  Transforms  (FFT)  and 
Periodograms)  are  typically  employed  to  determine  the  spectral 
content  of  the  data.  These  techniques  are  not  adequate  for 
time  series  data  samples  that  are  limited  in  number  (short). 

A  record  of  data  is  short  if  the  number  of  samples  from  the 
process  is  more  than  an  order  of  magnitude  smaller  than  the 
reciprical  of  the  lowest  frequency  of  interest.  For  example, 
a  short  period  recording  of  a  seismic  event  may  produce  less 
than  100  samples.  Seismological  events  often  occur  at  the 
lower  end  of  the  very  low  frequency  spectrum  (1  to  30,000 
Hertz).  Another  example  is  a  set  of  data  indicating  doppler 
frequency  shifts  based  on  a  very  few  (less  than  100)  radar 
returns.  This  thesis  considers  several  spectral  estimation 
techniques  and  determines  what  techniques  are  applicable  to 
analyze  short  record  time  series  data.  The  criteria  measur¬ 
ed  are  the  statistical  properties  and  computational  efficiency 
of  each  technique. 

The  spectral  estimation  techniques  considered  in  this 


thesis  are  Fast  Fourier  Transform,  smoothed  periodograms , 
Maximum  Entropy  Method  (MEM),  and  bandlimited  extrapolation. 
Methods  of  determining  the  statistical  properties,  i.e. 
mean  value,  autocorrelation  function,  and  variance,  of  the 
spectral  estimates  are  also  considered.  Additionally,  eval¬ 
uation  of  the  bias  in  the  statistical  estimates  is  pro¬ 
vided  along  with  a  determination  of  the  impact  of  the  bias 
on  the  statistical  estimates.  The  time  required  to  imple¬ 
ment  the  analysis  algorithms  (computer  execution  time)  and 
ease  of  application  are  also  considered  in  this  thesis. 

To  carryout  the  analysis,  several  basic  assumptions  are 
required.  Since  the  statistics  of  the  input  data  are  total¬ 
ly  unknown,  we  assume  that  the  data  are  samples  from  a 
stationary  random  process,  the  result  of  some  linear  system 
being  driven  by  an  unknown  noise.  Through  this  research, 
equal  spacud  samples  of  the  random  process  is  assumed. 

These  assumptions  are  justified  by  the  fact  that  many  of  the 
signals  of  interest  are  characterized  by  the  foregoing  as¬ 
sumptions  (Refs  2,7,24,25,28,  and  39).  The  intent  of  this 
research  is  to  provide  a  means  of  estimating  the  spectral 
density,  not  identification  of  the  spectral  density  data. 
Therefore,  hypothesis  testing  techniques  are  not  considered 
or  attempted. 

The  following  sections  provide  a  mix  of  current  theore¬ 
tical  knowledge,  general  application  techniques,  and  specific 
problem  analysis.  Section  II  lays  the  general  groundwork 
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in  spectral  estimation  techniques.  Methods  of  determining 
the  statistical  properties  of  the  estimated  spectrum  are  in¬ 
cluded.  Traditional  and  smoothed  periodogram  spectral  es¬ 
timation  techniques  are  discussed  in  this  section.  A  detail¬ 
ed  theoretical  background  of  MEM  is  included  along  with  three 
methods  of  determining  the  order  of  the  autoregressive  (AR) 
model.  Bandlimited  extrapolation  spectral  estimation  tech¬ 
niques  are  introduced.  In  particular,  a  method  developed  by 
Papoulis  using  the  FFT  is  suggested.  Section  II  concludes 
with  a  tabular  comparison  of  the  above  spectral  estimation 
techniques. 

In  Section  III  the  Bartlett  smoothed  periodogram  is  used 
as  a  "calibration  factor"  to  evaluate  the  averaged  results 
from  MEM  and  Papoulis  spectral  estimation.  Confidence  in 
these  techniques  is  derived  by  the  comparison.  Several 
known  input  signals,  including  simple  AR,  autoregressive  mov¬ 
ing  average  (ARMA),  and  sum  of  sinusoids  plus  Gaussian  noise, 
are  analyzed.  Finally,  each  input  process  is  evaluated  as  a 
short  record  (128  samples  or  less  with  a  sampling  interval 
of  one  second)  to  determine  analysis  results  under  near  real¬ 
istic  conditions. 

The  final  analysis  section  of  this  thesis  (Section  IV) 
provides  the  results  of  spectral  estimation  of  two  "real 
world"  problems.  The  MEM  and  Papoulis  techniques  are  imple¬ 
mented  to  determine  the  spectral  content  of  the  unknown  time 
sampled  signals.  A  comparison  of  the  resulting  data 
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indicates  that  the  MEM  technique  is  superior  to  the  constrain¬ 
ed  Papoulis  technique  when  frequency  resolution  is  a  prime 
factor. 

Section  V  concludes  the  research  indicating  the  utility 
of  the  MEM  and  Papoulis  spectral  estimation  techniques  for 
short  record  of  input  data.  Recommendations  for  future  re¬ 
search  in  this  area  are  suggested. 
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II .  Spectral  Estimation  Techniques 


Several  approaches  are  available  to  determine  the  power 
spectral  density  of  a  random  process.  The  methods  selected 
must  provide  accurate  results  when  only  a  limited  number  of 
samples  are  available.  Additionally,  the  complexity  of  the 
spectral  analysis  hardware  and  software  must  be  considered. 
This  section  outlines  the  techniques  considered  for  the  spec¬ 
tral  analysis  problem  and  provides  the  rationale  used  to 
select  the  maximum  entropy  spectral  estimation  method  and 
the  Papoulis  algorithm  as  the  techniques  for  the  problem 
analysis.  The  methods  initially  considered  are  as  follows: 

1.  Filter  banks  or  Fast  Fourier  Transforms  (FFT). 

2.  Smoothed  Periodograms  using  descrete  Fourier  Trans¬ 
forms  ( DFT ) . 

3.  Maximum  Entropy  Method  (MEM)  using  the  Burg  coef¬ 
ficient  estimation  techniques. 

4.  Papoulis  Bandlimited  Extrapolation  Spectral  Estim¬ 
ation. 

An  overview  of  these  techniques  is  provided  in  the  following 
subsections.  A  measure  of  accuracy  is  developed  through 
analysis  of  the  statistical  properties  of  the  estimated 
spectrum. 

Statistics  of  the  Estimated  Spectrum 

The  appropriate  techniques  to  estimate  the  spectral 
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density  of  an  unknown  input  process  are  developed  in  the  fol 
lowing  subsection.  The  accuracy  of  the  estimation  techni¬ 
ques  must  be  determined.  Just  how  good  is  the  estimate? 

Two  important  values  provide  an  indication  of  the  accuracy. 
The  mean  and  variance  of  the  estimated  spectral  density  in¬ 
dicates  the  true  value  of  the  spectrum  and  how  well  the  es¬ 
timate  approximates  this  value.  Two  measures  for  accuracy 
of  an  estimate  are  bias  and  consistency.  An  estimate  is 
biased  if  the  parameter  being  estimated  minus  the  expected 
value  (mean)  of  the  estimate  is  not  zero.  If  the  result  is 
zero  (unbiased),  the  estimated  value  is  the  true  value  (on 
the  average).  Therefore,  an  unbiased  estimate  is  more  ac¬ 
curate  than  a  biased  estimate.  An  estimator  is  said  to  be 
consistent  if  the  bias  and  the  variance  of  the  estimate 
tend  to  zero  as  the  number  of  observations  become  large 
(Refs  29:536  and  41:71).  In  other  words,  the  estimator  ap¬ 
proaches  the  true  value  as  the  number  of  observations  in¬ 
creases.  The  ideal  estimator  is  the  one  that  produces  es¬ 
timates  that  are  both  unbiased  and  consistent. 

The  estimated  mean  of  the  estimated  spectral  density 
is  calculated  using  the  sample  mean  equation  (Ref29: 536-537) 

N-l 

E(p(»n)>  -  v»n)  *  i?  I  ‘VV  O) 

i-0 

/\ 

where  Pi(un)  are  the  n  estimated  power  spectral  density 
values  for  the  ith  input  sequence  and  N  is  the  number  of 
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statistically  independent  sequences.  Assuming  that  mp(u)n) 
is  a  sequence  of  independent  Gaussian  random  variables,  the 

A 

density  function  of  mp(u)n)  is  also  Gaussian.  Therefore, 
the  expected  value  of  [®p(wn)J2  is  equal  to 


1 

W 


N-l  N-l 


I  E 

i«0  1=( 


E{Pi(o)n)Pj(a)n)} 
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and  the  variance  of  m  )  is 

P  D 


Var  mp(<*)n)j  ■  E{  mp(o)n)  2}-  E{mp((i)Q)  }j 2 


(3) 


From  Eqs  (2)  and  (3),  as  N  increases,  the  variance  of  the 
sample  mean  goes  to  zero  and  the  sample  mean  approximates 
the  actual  spectrum  arbitrarily  closely. 

The  maximum  likelihood  variance  estimate  of  the  estimat 
ed  spectrum  is  called  Var{p(wn)}  by  (Ref29:536-537) 


n-i  i  r  1 

Var{p(a)n)}=  J  £  VV  *“  V'V  * 
i-0L  J  J 


(4) 


and  is  shown  to  be  biased  by  Maybeck  (Ref27:398).  The  un¬ 
biased  variance  estimate  is  derived  and  is  shown,  for  the 
case  of  independent  observation,  by  Uaybeck  to  be 


*> 


N-l 

Var'(J(o.n))  -  jpj  £ 

i=0 

where  the  prime  (")  denotes  the  unbiased  variance  estimator. 
Thus,  if  N  is  large  enough,  N/(N-1)=1  ,  the  maximum  like¬ 

lihood  estimate  is  unbiased.  Also,  note  that  the  estimate  is 
consistent,  the  variance  estimate  tends  to  zero  as  N  is  in¬ 
creased.  Therefore,  if  independent  sequences  are  assumed, 
the  mean  and  variance  of  the  estimated  spectrum  is  determined, 
providing  information  concerning  the  bias  and  consistency  of 
the  estimate. 

The  following  paragraphs  provide  an  overview  of  the 
spectral  estimation  techniques  considered  in  this  research. 

Filter  Bank  or  Fast  Fourier  Transform 

A  straightforward  approach  to  spectral  analysis  is  to 
present  the  time  domain  signal  to  a  bank  of  narrowband  band¬ 
pass  filters.  Each  filter  will  only  pass  a  specific  range 
of  frequencies.  The  bank  of  filters  approximate  the  Fourier 
Transform  in  the  limit  as  the  number  of  filters  increases 
and  the  range  of  frequencies  decreases  (Ref  20:4). 

00 

x(u)  *  t)e"Ju,tdt  (6) 

•00 

The  technique  offers  almost  ideal  spectral  estimation  analysis, 
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limited  only  by  the  number  and  bandwidth  of  the  filters  and 
the  roll  off  characteristics  of  the  non  ideal  filters  actual¬ 
ly  implemented.  Typically,  the  system  is  realized  digitally 
and  is  usually  limited  by  computation  time.  The  amount  of 
computation  time  to  transform  a  given  set  of  time  data  to 
the  frequency  domain  is  inversly  proportional  to  the  band¬ 
width  of  each  filter  [the  Heisenberg  Uncertainty  Principle 
(Ref20:5)].  The  bandwidth  of  the  filters  determine  frequen¬ 
cy  resolution  (the  concentration  of  a  spectral  estimate  ex¬ 
pressed  in  frequency  units),  indicating  that  finer  resolu¬ 
tion  is  achieved  by  reducing  the  bandwidth  of  each  filter. 
Therefore,  if  the  input  is  finite,  and  short  as  described 
in  the  introduction,  frequency  resolution  may  be  jeopardiz¬ 
ed  when  computation  time  is  constrained. 

Another  inherent  difficulty  with  filter  bank  analysis 
is  the  window  function  effects.  The  finite  data  is  inter¬ 
preted  as  a  multiplication  of  the  original  signal  by  a  window 
function  (Ref  20:5).  The  result  is  a  spectrum  which  is  the 

convolution  of  the  original,  desired,  spectra  by  a  transform 
of  the  window  function.  Incorrect  selection  of  the  window 
function  may  result  in  analysis  that  is  truly  very  different 
than  the  actual  spectra. 

The  filter  bank  approach  is  thus  limited  by  the  trade¬ 
off  between  time  measured  data  and  frequency  resolution  ob- 
tainable,  and  the  proper  selection  of  the  window  function. 
Oppenheim  and  Schafer  (Ref 29: 549-553)  provide  an  excellent 
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discussion  of  the  effects  of  window  functions  and  the  rela¬ 
tion  to  correct  spectral  estimation.  Simply  stated,  win¬ 
dowing  is  the  result  of  convolving  the  data  of  interest  with 
some  window  function.  The  function  may  be  rectangular,  tri¬ 
angular,  and  other  shapes.  The  result  of  the  convolution, 
in  effect,  truncates  the  input  sequence  outside  the  period  of 
interest  and  conditions  the  data  inside  the  period.  The  win¬ 
dow  function  is  selected  to  insure  a  more  accurate  represen¬ 
tation  of  the  spectral  density  by  negating  the  effects  of 
the  FFT.  Close  relationships  exist  between  the  windowing  of 
the  original  input  data  and  periodogram  spectral  estimation 
techniques.  This  technique  is  discussed  in  the  following 
subsect ion. 

Smoothed  Periodograms  Using  Discrete  Fourier  Trans forms ( DFT ) . 

If  the  input  process  is  stationary  (initial  assumptions 
Section  I)  and  the  autocorrelation  function  is  known,  then 
the  power  spectral  density  of  the  input  can  be  determined 
(Ref30 : 338) . 

00 

P(u)  =  f  <Kt)e"Jwtdt  (7) 


where  the  Fourier  Transform  P(oj)  is  the  spectral  density 
and  <fr(t)  is  the  autocorrelation  function.  In  discrete  form, 

00 

p(«)  -  Ys  vk)e’Jwk  (8) 

k— 00 
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where  N  is  the  number  of  nonzero  input  autocorrelation  val¬ 
ues.  If  x(k)  is  a  real,  finite  length,  input  sequence 
0£k<N-l  then, 

N-l 

Vw)  =  £  x(k)e_ja)k  (9) 

k=0 

is  the  Fourier  Transform  of  the  input.  Since  the  autocorre¬ 
lation  function  of  the  input  sequence  is  unknown,  an  estimat¬ 
ed  value  is  introduced  (Ref29:542). 

N-| n| -1 

$n<k)  5  iRnJ  Z  |  x(k)x(k-n)  |  (10) 

k=0 

Combining  Eqs  (8), (9),  and  (10)  leads  to  an  estimate  of  the 
power  spectral  density. 


y<v  *  * 


w 


(id 


A 

where  Pp(<i>n)  denotes  a  power  spectral  estimate  using  Per- 
iodograms.  This  spectrum  estimate  is  called  a  periodogram. 

It  can  be  shown  that  the  periodogram  spectral  density 
is  biased  and  not  consistent  ( Ref 29: 540-545 ) .  Therefore,  al¬ 
ternative  techniques  must  be  considered  to  estimate  the  spec¬ 
tral  density.  Of  particular  importance  in  the  estimation 
technique  is  the  length  of  the  lag  for  each  estimate  with 
respect  to  the  total  sample  size  (Number  of  total  samples) 
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and  the  particular  window  function  used  in  the  frequency 
transformation.  Bartlett  and  Welch  (Refs  8  and  44)  consid¬ 
er  these  problems  in  the  development  of  smoothed  periodo- 
grams . 

One  technique  used  to  modify  the  periodogram  spectral 
estimator  so  that  the  estimate  is  consistent  is  attributed 
to  Bartlett  (Ref  8).  The  technique  is  particularly  use¬ 
ful  with  a  large  number  of  input  samples.  The  input  sequence 
is  segmented  into  several  sections  each  with  a  fixed  number 
of  samples.  The  spectral  estimated  from  all  sections  are 
averaged  to  produce  an  estimate  that  is  consistent.  The  es¬ 
timate  is  consistent  if  and  only  if  the  number  of  sections  is 
large.  If,  for  example,  the  input  consists  of  N  samples 
and  is  separated  into  K  sections  of  M  samples  each,  Eq(ll) 
is  applied  to  each  section  to  produce  K  periodograms.  Each 
periodogram  is  assumed  to  be  independent.  An  average  over 
all  the  periodograms  produces  an  estimate  of  the  spectrum. 

The  spectral  estimator  is  written  as: 

K 

V“n>  *  Z  I  5i<“n>  <12> 

i=l 

*  ^ 

where  P,(<u)  is  the  ith  periodogram  estimate  and  Pd(<d  )  is 
x  a  n 

the  Bartlett  smoothed  periodogram  spectral  estimate.  The 

A 

variance  of  P_(io  )  is  inversly  proportional  to  the  number 
d  n 

of  periodograms,  K  (Ref29:548).  Therefore,  as  K  becomes 
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large,  the  variance  approaches  zero  and  the  Bartlett  esti¬ 
mate  is  consistent.  Additionally,  the  bias  of  the  Bartlett 
estimator  is  greater  than  that  of  the  simple  periodogram 
(Ref  8).  These  two  results  indicate  that  the  reduction  of 
the  variance  is  accomplished  at  the  expense  of  increased 
bias.  However,  increasing  K  results  in  fewer  samples  for 
each  periodogram,  thereby  decreasing  frequency  resolution. 

For  example,  if  it  is  known  that  the  spectrum  has  a  narrow 
peak,  M  must  be  chosen  large  enough  to  resolve  the  peak. 

For  small  variance,  K  must  be  large.  The  number  of  sam¬ 
ples  of  the  input  sequence  is  determined  by  the  variables 
K  and  M  .  Clearly,  if  low  variance  and  fine  frequency 
resolution  are  implicit,  an  extremely  long  input  record  may 
be  required. 

Blackman  and  Turkey  suggest  another  approach  to  spectral 
estimation  using  smoothed  periodograms  (Ref  9).  The  techni¬ 
que  is  to  smooth  the  periodogram  by  convolution  with  an  ap¬ 
propriate  spectral  window.  Welch  applied  the  Blackman-Tur- 
key  hypothesis  to  the  averaged  periodogram  solution  and  de¬ 
veloped  the  spectrum  estimation  technique  of  averaging  mod¬ 
ified  periodograms  (Ref  44).  Welch  modifies  the  Bartlett 
procedure  by  applying  a  window  function  directly  to  the  data 
input  sections.  The  periodogram  (Eq(ll))  is  thus  modofied 

to„  « 

*pw<“n>  ‘  Y. 

n«0 
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where 


x^*\n)  is  the  nth  input  sample  from  the  ith  section. 


and 


M-l 

D  "s  Z  *2<n>  <14> 

n=0 


where  W(n)  is  the  desired  window  and  P  ( u>  )  is  the  win- 

pW  D 

dowed  periodogram  function.  The  smoothed  spectral  estimator 
is 

K 

P  (w  )  =  i  V  p  (0)  )  (15) 

wv  n'  K  /_  pwv  n'  v  ' 

i=l 

Welch  shows  that  with  the  inclusion  of  U  ,  the  estimate 

A 

Pw(u>n)  is  not  biased  and,  if  the  sections  do  not  overlap, 
the  variance  is  inversly  proportional  to  K  .  Therefore,  the 
Welch  estimator  is  consistent.  Clearly,  the  Bartlett  pro¬ 
cedure  is  a  simplification  (subset)  of  Welch's  method  with 
the  window  function  W(n)  equal  to  one. 

For  the  specific  application  under  consideration  in 
this  thesis,  smoothed  periodograms  are  not  practical.  For 
accurate  spectral  estimation,  these  procedures  must  have  a 
large  data  record.  Bartlett  smoothed  periodogram  techniques 
do  provide  a  means  of  verifying  the  results  of  analysis  on 
known  signals  under  known  conditions  when  the  MEM  and  Pap- 
oulis  techniques  are  verified  in  Section  III. 
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Maximum  Entropy  Spectral  Estimation 


The  concept  of  linear  prediction  is  based  on  the  assump¬ 
tion  that  any  linear  system  can  be  represented  in  discrete¬ 
time  by  (Ref25:562) 


M  P 

x(n)  =  -  ^  0kx(n-k)+A  ^  Yjy(n-j ) , y0=l 
k=l  j=0 


(16) 


where  y(.)  is  some  unknown  input,  8k  ,  Yj  ,  and  A  are 
some  unknown  parameters,  and  l^k^M  and  lijip  .  In  other 


words,  the  output,  x(n)  ,  is  predicted  by  some  linear  com¬ 
bination  of  previous  outputs  and  inputs.  A  special  case  of 
this  linear  model  is  the  autoregressive  (AR)  (all-pole)  model 
defined  by 


MI 

x(n)  =  3kx(n-k)+y(n) 


(17) 


and  M  is  the  system  order.  This  equation  is  the  model  as¬ 
sumed  to  be  valid  for  the  analysis  under  consideration.  The 
equation  can  be  interpreted  as  a  stochastic  difference  equa¬ 
tion  and  describes  the  samples  of  a  sequence  from  a  random 
process  where  y(n)  is  a  white  noise  sequence  with  zero 
mean  and  variance  a*  .  Finite  order  AR  processes  are  use¬ 
ful  in  representing  many  physically  significant  linear  ran¬ 
dom  processes  (Ref  1,11,12,24,  and  25).  J.P.  Burg  provides 


the  power  spectral  equation  for  this  AR  Uodel  (Ref  12). 


$(u>m) 

m 


PM+lAt 


M  2 

-ju>mnAt 


1*.  I  Sne 

n=l 


(18) 


where  At  is  the  input  sampling  interval.  The  proper  values 
for  the  prediction  coefficients  ,  the  prediction  error 

power  PM+1  >  and  the  model  order  M  must  be  determined. 
These  values  are  discussed  in  the  following  paragraphs. 

The  relationship  between  the  AR  process  and  the  maxi¬ 
mum  Entropy  Method  is  extensively  researched  (Refs  13,38,39 
and  40).  A  brief  overview  of  the  relationship  is  provided 
here. 

Consider  a  situation  where  M  different  things,  , 

happen  with  probability  P.^  .  Information  is  defined  as 

-log^j^  where  P.^  is  a  particular  event  (Ref  47:428).  If 
all  the  Pi  are  known  to  be  equal  a  priori,  no  Information 
is  available.  However,  when  a  particular  PA  is  known,  a 
certain  amount  of  Information  is  gained.  Clearly,  the  prob¬ 
ability  of  occurence  is  related  to  Information.  Shannon  (Ref 
36)  explored  this  Information  Theory  to  measure  the  average 
Information  required  to  transmit  a  given  message.  Shannon 
calls  the  measure  entropy  and  is  the  average  Information  per 
time  interval. 

M 

H  ■  -k  V  P.  loe  P.  (19 


where  k  is  a  constant  and  logr  is  a  logarithm  base  r. 
Consider  a  process  xt  that  can  take  on  the  values  xi,X2,... 
xQ  .  Assume  the  available  Information  about  xt  is  the  aver 
age  values  <fi(xt)>  ,  <f2(xt)>  <fm(xt)>  of  some 

functions  fi(xt),  where  m-cn  .  The  probability  distribu¬ 
tion  Pt  =  F(xt)  that  is  consistent  with  the  Information 
and  is  maximally  free  of  other  limitations  maximizes  the  en¬ 
tropy  Eq( 19) .  In  the  continuum  case  (Ref  38:413), 


( 


H 


/I(2 N>  ln(of(£N>}d5 


(20) 


where  f(zN)  is  a  joint  probability  density  for  N  elements, 
and  C  is  a  constant. 

Smylie  et  al  (Ref  38:410-414)  show  that  the  entropy 
measure  H  is  not  a  valid  measure  for  frequencies  and  sug¬ 
gest  an  entropy  rate  h  that  relates  the  power  spectral  den¬ 
sity  to  entropy  rate.  The  suggested  entropy  rate  for  a 
stationary  Gaussian  process  is  defined  in  terms  of  the  spec¬ 
tral  density  (Ref38:414) 


h 


In  P(u))du> 


(21) 


where  f^  is  the  Nyquist  frequency  and  In  is  the  natural 
logarithm.  Replacing  P(<u)  with  the  power  spectral  density, 
Eq  (8),  yields 
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(22) 


assuming  a  uniform  sampling  rate  and  sampling  interval  of  one. 
Smylie  et  al  (Ref  38:414)  continued  the  derivation  to  yield 

h  =  i  In  (det  [*J  }  (23) 

wb^re  is  the  positive  semidefinite  M  by  M  matrix  of 

autocorrelation  coefficients.  Smylie  shows  that  4>„  is  equi- 
diagonal  and  indicates  that  such  matricies  are  named  Toeplitz. 

To  maximize  the  entropy,  must  be  determined  in  such 

a  way  as  to  impart  no  assumptions  about  the  input  data.  No 
erroneous  information  is  added  to  the  analysis.  This  restraint 
departs  extensively  from  the  previous  estimation  techniques. 

In  the  Bartlett  technique,  the  data  outside  the  window  is  as¬ 
sumed  to  be  zero,  suggesting  that  the  autocorrelation  function 
is  zero  outside  some  limiting  time  period.  If  the  autocorre¬ 
lation  is  zero,  periodic  extensions  used  with  FFTs  are  im¬ 
plicit.  Consequently,  additional  data  is  erroneously  added. 

As  a  result,  the  estimated  spectral  data  may  be  somewhat  dif¬ 
ferent  than  the  actual  spectrum.  The  Maximum  Entropy  tech¬ 
nique  provides  an  alternate  choice  for  selecting  autocorre¬ 
lations,  one  in  which  no  Information  is  added  or  implied. 

The  MEM  technique  implies  is  calculated  in  such  a  way 

as  to  maximize  the  uncertainty  or  entropy  and  make  no 
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assumptions  about  the  data  outside  the  time  sample 


To  calculate  the  MEM  power  spectrum[Eq(18)]  ,  the  order 
of  the  AR  process  (length  of  the  prediction  filter)  and  the 
prediction  coefficients  (BQ)  must  be  determined.  Procedures 
to  determine  these  values  are  discussed  in  the  following  para¬ 
graphs  . 

The  toeplitz  autocorrelation  matrix  for  an  Mth  order  pro¬ 
cess  is 


<J>(0)  <J>(1)  .  .  .  <KM) 

4>(1)  <t>(0)  .  .  .  <f>( M-l) 


(24) 


4>(M)  <KM-1)  .  .  <K0) 


if  it  is  assumed  that  4>(0)  ,<K1)  , . . .  ,<KM)  are  known.  Sup¬ 
pose  4>M+1  is  to  156  estimated  from  this  data,  then  <J>(M+1) 
must  be  selected  to  maximize  the  entropy.  Also.  <|>QI+2)  and 
so  on.  Therefore,  from  Eq(23),  <|>(M+1)  is  determined  by  max 

imizing  det  [<1>M+1]  with  respect  to  <KN  +1)  .  Since  $M+1 

must  be  positive  semidef inite,  det[$M+j]  has  a  single  max¬ 
imum  when  the  product  rule  of  differentiation  is  applied. 

With  4>m+1  determined,  <f>M+2  can  be  determined  in  a  similar 
manner.  For  the  case  maximizing  det^fl)  M+1J  with  respect  to 
$(M+1)  the  result  is 


4>(M+l)  4>(M)  .  .  .  <K1) 


Solving  Eq(17)  for  AR  order  M  , 

x(n)=eix(n-l)+B2x(n-2)+  ...  +3Mx(n-M)+y(n)  (26) 

Multiplying  through  by  x(n-k),k«0,  and  taking  expected  values, 

4>(k)  =  3 1  (j>(k-l)  +  B2  <()(k-2 )+  ...  +g^(k-M)  (27) 

where  E{x(n-k)y(n) }  =  0  for  a  zero  mean  Gaussian  input  y(n). 
Substituting  the  values  of  k=l,2, . . . ,M+1  into  Eq  (27) 
yields  a  set  of  simultaneous  equations. 

4.(1)  -3 14>(0)-  .  .  .  -3m4>(M-1)  =  0 

•  •  •  • 

(28) 


<KM+1)  -  3 i 4>(M)-  .  .  .  -3m4»(1)  =  0 

These  equations  are  known  as  the  Yule-Walker  equations  (Refs 

25, 43, and  46).  If  4>(0) , . . . , 4»(M)  are  known,  3i . 3M 

can  be  determined.  Also,  4>(M+J.)  can  be  determined  by 


solving  for  the  determinant  of  the  matrix  of  all  the  <f>(i) 
set  to  zero. 


♦(1)  4(0)  • .  .  .  4(M-1) 

4(2)  <(>(1)  .  .  .  4 ( M-2 ) 


=  0 


(29) 


(KM+1)  <KM)  .  .  .  4(1) 


but  this  equation  is  identical  to  Eq(25).  So  the  Yule-Walker 
equations  are  identical  to  maximizing  the  entropy  of  the  pro¬ 
cess. 

Suppose  k  [Eq(27)J  is  allowed  to  equal  zero.  The  ex¬ 
pected  value  of  x(n-k)y(n )  is  E{x(n-k)y(n) }=  E{x(n)y(n)}  , 
but  E{x(n)y(n)}=  E{y2(n)}*a2y  for  the  zero  mean  Gaussian 
model  Eq  (17).  Eq  (27)  is  now,  for  k=0, 

<K0)-e1<Kl)-...-eM<J>(M)«a2  (30) 

and  Eq  (28)  can  be  augmented  to  include  all  k=0,l,...,M  as, 


■<j>(0)  4(1)  ....  <j,(M) 

1 

r  a2  ■ 

<K  1)  4.(0)  ....  <t>(M-l) 

•  •  • 

• 

•  •  • 

• 

• 

• 

y 

0 

• 

• 

(31) 

•  •  • 

• 

4(M)  <{>(M-1)  .  .  .  4(0) 

• 

--V 

• 

0 
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which  is  an  alternate  form  of  the  Yule-Walker  equations.  In 
general,  the  variance,  a*  ,  is  estimated  as  the  prediction 

O' 

error  power  PM+1  ,  resulting  from  the  convolution  of  x(n) 
with  the  M+l  point  prediction  error  filter  Eq(17)  (Ref 39: 
187).  is  often  referred  to  as  the  innovations  power. 

If  the  data  are  known  for  the  M-l  case,  the  data  for 
the  Mth  case  is  determined  by  calculating  the  coefficients 
6i,...,Bu,  autocorrelation  4> ( M) ,  innovations  power  Pu+1, 
and  then  calculating  the  power  spectrum  Eq(18).  For  the 
M+l  case,  there  are  M+l  equations  and  M+2  unknowns. 

Burg  (Refs  11,12  and  13)  develops  an  elegant  procedure  that 
does  not  make  any  assumptions  about  the  data  outside  the  sam¬ 
ple.  Additionally,  the  technique  does  not  require  a  prior 
estimate  of  the  autocovariance  or  autocorrelation  function. 
Burg  suggests  that  the  error  power  be  calculated  by  running 
a  prediction  error  filter  over  the  data  in  a  forward  and 
backward  direction,  but  not  off  the  data.  As  each  set  of  co¬ 
efficients  is  estimated,  the  error  filter  order  is  increased 
by  one  and  the  procedure  is  repeated  until  the  desired  model 
order  is  reached.  The  appropriate  model  order  is  discussed 
in  later  paragraphs.  The  recursion  corresponds  to  an  approx¬ 
imate  maximum  likelihood  estimation  of  the  AR  coefficients 
(Refs  10:46-84  and  13:375). 

Consider  the  case  of  a  set  of  N  input  samples  of  data, 
x(l),  ...,  x(N),  all  with  equal  spacing  At  ,  and  the  power 
spectrum  of  an  Mth  order  AR  model  is  required.  First, 
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assume  the  model  order  is  zero.  Therefore,  an  estimate  of 
the  autocorrelation  for  lag  one  is 


M1)  R  Y  x2(n) 


(32) 


where  the  (1)  in  $0(1)  denotes  the  first  lag.  From  Eq(31), 
Pi*4>o(0)  .  The  system  for  model  order  one  is 


1  4>o(l) 

$o(D  Pi 


(33) 


where  Bn  is  the  estimate  of  Bi  .  The  forward  prediction 


error  is,  from  Eq(17) 


ei (n)  *  x(n)  -  x(n) 


x(n)-Bi ix(n-l)  ;n=2.3. . . . ,N(34) 


Similarly,  the  backword  prediction  error  is 


ei(n)  *  x(n)-Bnx(n+l)  ;  n=N-l,N-2, . . . ,1 


(35) 


The  prediction  error  power  is  calculated  by 


Pl  “  STn-1')  Y.  {L*(n+l)-Bnx(n)] 


+[x(n)-Bi ix(n+l)J  2 ) 


(36) 
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where  Bn  is  chosen  to  minimize  Pi  .  Bn  is  calculated 


by  (Refs  5:70  and  11:5) 


*1 

2  ^  x(n)x(n+l) 


n=l 

N-l 


(37) 


^  [x2 (n)+x2 (n+1)] 


Solving  for  $(1)  and  P2  yields 


i  ( l )  =  6i  i  d>  ( o ) 


(38) 


P2  =  Pi(l-6n2) 


(39) 


In  general,  for  the  step  from  model  order  M-l  to  M  , 

(40) 

N-M  L  M  2  M  i 

PM  ~  S(M-10  Z  |x(n)-^  SMJx(n-j)l  +  bt(n+M)- ^  BMJx(n+M-J )] 
n=i  L  j=l  J  *■  j=i  J 

where  PM  is  minimized  with  respect  to  the  estimated  coeffi¬ 
cient  6mm  (Ref5:71).  Calculation  of  the  prediction  coeffi¬ 
cients  ewk  ;  j=l ,2 , . . . ,M-1  is  accomplished  by  the  Levinson- 
Durbin  procedure.  Makhoul  (Ref25:566)  outlines  the  procedure 
which  is  a  technique  for  solving  for  the  gMj  data  with  esti¬ 
mated  $(j)  and  Eq(31).  The  procedure  solves  for  the  data 
in  only  M2  operations  and  requires  2M  storage  locations. 
The  recursion  equation  for  fju.  is 

MJ 
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'Mj  M-lj  MM  M-lM-j  '■ 


From  Eq(31), 
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where  x  indicates  a  don’t  care  condition,  and 


Finally, 


M 

*m(m>  *  -  <44> 

j-1 

Anderson  (Ref  5)  provides  an  algorithm  to  carryout  the  recur¬ 
sions.  The  routine  is  suitable  for  implementation  on  a  com¬ 
puter  and  is  used  in  the  analysis  sections  of  this  report. 

The  coefficients  3  jji *  •  •  •  and  the  err°r  power  (innova¬ 

tions  power)  are  implemented  in  Eq(18)  to  produce  an  estimate 
of  the  power  spectrum  of  AR  model  order  equal  to  M 

Statistical  properties  of  the  MEM  spectral  estimate  are 
investigated  by  Makhoul  and  Kromer  (Refs  23  and  25:568-569). 

It  is  determined  that  the  MEM  estimate  is: 

1.  Asymptotically  unbiased,  as 
M-*-<®, 

E{P(u>n)}  =  P(u)  (45) 

2.  Asymptotically  Normal  with  variance 

Var{P(u,n)}  =  |?!{p2(a,n)}  (46) 

at  each  frequency 

Thus,  for  large  M  and  N  and  a  smooth  true  spectrum,  any 
spectrum  can  be  approximated  arbitrarily  closely  by  the  AR 
model  (Refs  39:192  and  25:569).  Fuller  makes  this  statement 
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even  stronger  in  a  theorem  for  any  continuous  spectral  den¬ 
sity  (Ref 18: 147-149). 

To  estimate  the  power  spectrum  using  MEM  ,  the  order 
of  the  AR  model  must  be  correctly  determined.  Three  techni¬ 
ques  are  suggested  in  the  following  paragraphs.  Selection 
of  the  technique  to  be  used  in  analyzing  unknown  data  in¬ 
puts  is  primarily  contingent  on  the  simplicity  of  implemen¬ 
tation. 

One  order  selection  method  considers  the  well-known 
statistical  analysis  tool,  the  Kolmogorov-Smirnov  single  sam¬ 
ple  test  for  goodness  of  fit.  The  test  is  found  in  several 
sources  (Refs  10,15,  and  37).  Additionally,  computer  pro¬ 
grams  are  generally  available  in  program  libraries. 

An  estimate  of  the  input  sequence  is  generated  by  apply¬ 
ing  the  prediction  coefficients  Eq(41) 

M 

x(n)  =  £eMjx(n-j)  (47) 

j  =  l 

for  model  order  M  ,  and  the  error  is  determined 

•  e(n)  =  x(n)-x(n)  ;  l£n£N  (48) 

l 

If  the  prediction  coefficients  are  correctly  estimated  and 
the  order  is  correct,  then  e(n)  equals  y(n)  [Eq(17)J. 

The  initial  assumption  is  that  the  system  [Eq(17)]  is  driven 
by  white-Gaussian  noise  ( WGN ) .  The  power  spectral  density  of 
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WGN  is  constant  over  all  frequencies.  Therefore,  E(m)  ,  the 
power  spectral  density  of  e(n)  ,  will  be  constant  over  some 

A 

bandwidth  if  x(n)  is  a  correct  estimate.  Integrating  E(w) 
over  this  bandwidth  will  produce  a  ramp  function  approximating 
the  uniform  distribution.  The  Kolmogorov-Smirnov  test  compares 
the  integrated  E( )  with  an  actual  uniform  distribution  func¬ 
tion,  and  generates  an  output  that  represents  a  probability  es¬ 
timate  of  how  well  the  integrated  E(ui)  approximates  the  uni¬ 
form  distribution.  This  probability  output  for  order  equal  to 
M  is  compared  to  previous  values  from  similar  test  on  orders 
up  to  M  .  The  maximum  Kolmogorov-Smirnov  probability  corre¬ 
sponds  to  the  best  estimate  of  the  system  order. 

The  Kolmogorov-Smirnov  test  clearly  necessitates  consider¬ 
able  manipulation  of  the  data.  Therefore,  the  test  requires 
a  large  amount  of  programming  and  computer  time.  Because  of 
these  factors,  Kolmogorov-Smirnov  tests  are  not  included  in 
the  final  analysis  section  of  this  paper.  However,  because 
of  the  confidence  many  statistical  practitioners  have  in  this 
test,  it  is  used  as  a  verification  and  confidence  building 
tool  for  the  following  two  model  order  determining  techni¬ 
ques. 

The  problem  of  model  order  determination  has  been  ex¬ 
tensively  investigated  in  the  past  few  decades.  Many  of  the 
techniques  result  in  considerable  computations  and  many  as¬ 
sumptions  (Refs  6,16,17,34,  and  45).  Two  techniques  develop¬ 
ed  by  Akaike  (Refs  1,2, 3, and  4)  show  considerable  promise 


28 


and  are  used  as  analysis  tools  to  evaluate  the  problem  con¬ 
sidered  in  this  paper. 

Akaike's  first  technique  adopts  a  decision  theoretical 
approach  in  which  a  value  of  merit  is  defined  for  each  AR 
model.  The  figure  of  merit  is  the  mean-square  error  of  the 
one-step-ahead  prediction  determined  by  using  the  MEM  to  es¬ 
timated  the  coefficients,  £  .  In  other  words,  if  x(n) 

is  defined  as  the  estimated  prediction  of  x(n)  ,  then 

FPE  =  E{  [x(n )-x(n)J  2 }  (49) 

where  FPE  is  called  the  final  predictor  error.  Another  esti¬ 
mate  (y(n))  is  generated  from  the  same  coefficients  and  an  in¬ 
dependent  realization  y(n)  .  FPE  for  this  case  is 

FPE  =  E{ [y(n )-y(n )]  }  (50) 

Akaike  (Ref  1).  shows  that  the  expression  for  the  FPE  of  the 
Mth  order  AR  process  is  then 


(FPE)m  = 


N+M+l 

N-M-l  fM+1 


(51) 


where  N  is  the  total  sample  size  and  PM+^  is  as  defined 
by  Eq(31).  The  scheme  is  to  determine  some  maximum  order 
L<N  .  Then  calculate  FPE,  Eq(51),  for  successive  order 
M(M=1,2 , . . . ,L) .  The  minimum  FPE  yields  the  estimate  of  the 
AR  model  order  for  preduct  ion.  FPE  order  estimation  techni¬ 
ques  are  investigated  by  several  authors  and  provide 
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excellent  results  (Refs  3, 22,39, and  40).  This  technique  is 
simply  implemented  and  requires  only  one  subjective  assump¬ 
tion,  namely  selection  of  the  maximum  possible  order  L 

In  a  recent  paper  by  R.H.  Jones  (Ref  22),  it  is  deter¬ 
mined  that  the  FPE  scheme  outlined  by  Akaike  becomes  biased 
as  the  model  order  approaches  the  sample  size.  This  result 
is  verified  by  experimental  investigation.  Jones  shows  that 
if  the  input  data  is  white  noise,  the  estimated  one  step  pre¬ 
diction  variance  remains  constant  as  model  order  increases 
to  N/2(for  an  N-sample  input)  .  For  orders  above  N/2  , 

the  variance  decreases  linearly  as  model  order  increases. 

With  white  noise  as  an  input,  the  variance  should  remain 
constant  for  all  models.  To  off-set  the  decreased  variance, 
Jones  suggests  the  following  modifications  to  the  FPE  Eq(51) 


(FPE)m 


N+M+l 

N-M-l 


N+M+l 

N-M-l 


PM+i  0<M<N/2 

Pm+i/U.S-M/N)  N/2<M<N-1 


(52) 


Applying  this  new  FPE  equation  yields  a  constant  variance  for 
all  model  orders  up  to  N-l 

In  subsequent  work  Akaike  extended  the  FPE  scheme  to  be 
applicable  to  any  maximum  likelihood  situation  (Ref  1).  He 
describes  an  information  criterion  for  estimating  the  model 
order  that  employs  the  ln-maximum  likelihood  estimate  of  a 
given  process.  Akaike’ s  A-information  criterion  (AIC)  is 
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defined  by 


(AIC)m  *=  -2  In  (maximum  likelihood)  +  2U  (53) 

"A"  denotes  the  first  case,  with  the  expectation  that 

may  follow.  It  turns  out  that  for  the  case  of  univariate  AR 

model  with  Gaussian  input  (Ref  21), 

(AIC)m  -  N  In  (Pm+1)+2(1+M)  (54) 

Again,  a  maximum  order  L<N  is  selected  and  the  minimum  AIC 
denotes  the  proper  estimate  of  the  model  order.  R.H.  Jones 
(Ref  21:895)  shows  that  AIC  and  FPE  select  the  same  order. 
Additionally,  Jones  suggests  that  FPE  and  AIC  are  asymptotic¬ 
ally  equivelent.  AIC,  like  FPE,  is  easily  implemented  and  is 
therefore,  a  desirable  analysis  tool  for  the  problem  under 
consideration. 

The  Maximum  Entropy  spectral  estimation  techniques  out¬ 
lined  in  the  preceding  paragraphs  provide  an  excellent  tool 
for  the  research  in  this  paper.  Since  the  problems  are  de¬ 
scribed  for  short  record  data,  MEM  is  ideally  suited.  Addi¬ 
tionally,  MEM  provides  an  unbiased  estimate  of  the  spectral 
density  of  the  input  sequence. 

Papoulis  Bandlimited  Extrapolation  Spectral  Estimation 

A  new  algorithm,  developed  by  A.  Papoulis  (Refs  31,  32, 
and  33)  appears  to  be  applicable  to  the  problem  under  con¬ 
sideration.  The  technique  involves  an  iterative  routine 
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using  only  a  discrete  FFT/IFFT.  Specifically,  Papoulis' 
algorithm  is  applicable  to  stationary  bandlimited  signals. 
The  algorithm  is  based  on  one  simple  assumption 

P(w)  ■  0  whenever  |w|>a  (55) 


Note  that  this  assumption  does  not  restrict  the  time  function 
x(n)  to  some  specific  value  outside  the  period  of  interest. 

The  problem  is  to  find  the  Fourier  Transform  of  a  band- 
limited  function  x(t)  in  terms  of  some  simple  segment  of 
x(t),  x(n)  .  Papoulis  suggests  that  the  solution  is  found 
by  the  following  iteration: 


1 .  Form  x i ( t ) 


x(t)  | t | <T 
0  | t I >T 


(56) 


2.  Determine  the  FFT 


xi(“)  “  £  [xi(t^ 


(57) 


3 .  Form 


P,(0)) 


(to)  |  (D  |  <0 
0  |  cu  |  >0 


4.  Compute 

x2(t) 


x'(t) 


'{P^w)} 
x(t)  | t | <T 

x2(t)  | t | >T 


(58) 


(59) 


32 


5.  Continue  steps  2,3,' and  4  to  the  nth  iteration 


to  yield 
pn(“>  = 


|  ui  |  <a 
|  oo  |  >a 


(60 


The  number  of  iterations  ,n,  can  be  controlled  by  a 
preset  mean— square  error  (MSE)  value.  The  MSE  for  the  nth 
iteration  is 


oo  O 

in  -  J[x(t)-xn(t)]2dt+->  ^  J  [P(«o)-Pn(M)J 


:  doj 


(63 


-o 


/ 


•  /' 


[x(t)-xn(t)]2dt  +  I  [x(t)-xn(t)]2dt 


t  <T 


t  >T 


rewriting  the  far  right  integral  in  terms  of  the  Fourier 
Transform 


J [X(t)-xn(t)]2  +  ||P(.)-Xn+1(U>l2 

1 1 1  <T  "]w|>o 


dco 


P(u))-Pn+1(u)  I  2du 


(6 


From  Eqs(61)  and  (62), 

en"en+l  =  f[x<t)-xn(t)]2dt+i¥  f  lP(w)'Xn+l(lo)  1 2du)  (€ 
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The  two  integrals  on  the  right  side  of  Eq(63)  are  never  less 

than  zero.  Therefore,  e  -e  ...  is  never  less  than  zero.  It 

n  n+l 

follows  then,  that  en  is  monotone  non-increasing.  This  con¬ 
clusion  indicates  that  the  iteration  procedure  can  continue 
to  some  value  to  satisfy  a  preset  condition. 

Papoulis  (Ref  31:738-740)  proves  the  convergence  of  the 
technique,  eQ-*-0  as  n-*-«  in  the  deterministic  case.  The 
proof  is  through  a  comparison  of  the  technique  to  a  method 
for  extrapolation  of  band-limited  functions  based  on  the  ex¬ 
pansion  of  xi(t)  into  a  series  of  prolate  spheroidal  func¬ 
tions.  The  results  are  provided  in  Appendix  A.  The  algorithm 
is  used  by  Gerchberg  (Ref  19),  in  a  1974  report  on  image  pro¬ 
cessing,  with  satisfactory  results. 


Comparison  of  Spectral  Estimation  Techniques 


\  i 


To  complete  the  background  information  concerning  tech¬ 
niques  for  spectral  estimation,  several  areas  must  be  consid¬ 
ered.  In  particular,  the  techniques  must  be  applicable  to 
analysis  of  short  input  record.  However,  some  aspects  of  the 
traditional  techniques  can  be  used  to  verify  the  newer  or  un- 
familier  techniques. 

Articles  by  Burnard  and  LaCoss  (Refs  7  and  24)  provide 
additional  data  and  solidify  the  conclusions  concerning  the 
various  techniques.  For  example,  Burnard  and  LaCoss  suggest 
that  the  variance  of  the  spectral  estimate  for  both  MEM  and 
Periodograms  is  proportional  to  2L/N,  where  L  is  the  number 


of  sections  from  the  total  of  N  samples.  Also,  LaCoss  shows 
that  frequency  resolution  is  inversely  proportional  to  N 
for  Periodograms  and  inversely  proportional  to  N2  for  MEM. 
These  results,  along  with  the  conclusions  derived  in  this 
section  are  provided  in  Table  I.  Clearly,  the  MEM  and  Pap- 
oulis  techniques  are  applicable  to  the  problem  under  consider- 
at ion . 

Several  different  input  sequences  are  applied  to  both 
MEM  and  Papoulis  spectral  estimation  computer  coded  algorithms. 
It  is  determined  that  on  the  average,  MEM  routines  complete 
execution  in  about  the  same  time  as  the  Papoulis'  routine 
when  the  number  of  iterations  is  less  than  100.  Therefore, 
the  number  of  Papoulis  iterations  is  limited  to  a  maximum  of 
100  to  allow  a  rough  equivalence  in  computational  time  be¬ 
tween  the  MEM  and  Papoulis  techniques. 

Smoothed  periodogram  techniques  provide  the  basis  of 
analysis  to  verify  the  MEM  and  Papoulis  techniques  applied 
to  known  signal  samples.  The  following  section  provides  the 
results  of  applying  the  techniques  derived  in  this  section. 
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Comparison  oi  Spectral  Estimation  Techniques 


o 
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Analysis  of  Known  Processes 


The  purpose  of  this  section  is  to  provide  graphic  illus¬ 
trations  of  the  spectral  analysis  techniques  developed  in  the 
previous  section.  Known  time  sampled  data  is  applied  to  the 
estimation  equations  to  determine  the  spectral  content.  The 
Bartlett  smoothed  periodogram  technique  provides  true  spectral 
density  data  for  verification  of  the  estimates  developed  by 
the  MEM  and  Papoulis  techniques.  Calculation  of  the  variance 
of  the  estimated  spectrum  provides  additional  information  a- 
bout  accuracy.  Three  specific  examples  are  evaluated  to  de¬ 
velop  confidence  in  the  MEM  and  Papoulis  techniques. 

The  techniques  developed  in  Section  II  are  applied  in  sev¬ 
eral  subroutines.  These  routines  have  been  coded  using  For¬ 
tran  IV  and  are  provided  in  Appendix  B.  Figure  1  shows  a  sim¬ 
plified  flow  diagram  of  an  analysis  routine  to  determine  the 
spectral  density  and  statistics  of  an  input  process.  The  se¬ 
quence  is  assumed  to  be  very  large  (greater  than  30,000  points), 
to  facilitate  accurate  determination  of  the  Bartlett  periodo¬ 
gram,  the  sample  mean,  and  the  estimated  variance  of  the  MEM 
and  Papoulis  techniques.  These  values  are  calculated  from 
64  independent  realizations  each  with  128  samples. 

Figure  2  illustrates  the  algorithm  used  to  determine  the 
estimated  spectral  density  using  MEM  and  Papoulis  techniques. 
This  routine  assumes  that  the  input  sequence  is  short  as 
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Initialize  Data  and  Plotter 
Plot  128  Samples  of  Input  Data 


Determine  FPE,AIC,  and  Kolmogorov- 
Smirnov.  AR  model  order 
Calculate  MEM  Coefficients 
.CALL  MEMC01 
.CALL  AFPE 
.CALL  AAIC 
.CALL  KOLSMR 

Loop  of  Orders  0-31  _  _ 


Plot  FPE,  AIC,  and  Kolmogorov- 
Smirnov  data  for  Orders  0-31 


z 


Decide  on  AR  Model  Order  / 

I 


Determine  Sample  Mean  and 
Variance  for  MEM  Estimate,  64  Windows 
(Skip  50  Data  Points  Between  Each  Window) 
Let  N=64  (64  Windows) 

.CALL  BTLET 


Plot  Bartlett  Spectrum 


Figure  1.  Flow  Diagram  to  Determine  AR  Model  Order,  Bartlett 
Spectrum,  Sample  Mean,  and  Variance  of  MEM  and 
Papoulis  Estimation 
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Plot  Bartlett  Spectrum  / 


1 

Determine  Sample  Mean  and 
Variance  for  MEM  Estimate,  64 
Windows  (Skip  50  Data  Points 
Between  Each  Window) 

(Eqs(l)  and(4)) 

.CALL  MEMC01 
.CALL  PSPEC 

Loop  64  Times _ 


Plot  Sample  Mean  and  Variance^ 


Determine  Sample  Mean  and 
Variance  for  Papoulis  Technique, 
64  Windows  (Skip  50  Data  Points 
Bewteen  Each  Window) 

(Eqs(l)  and  (4)) 

.CALL  PAPO 


Loop  64  Times 


Plot  Sample  Mear^and  Variance  / 


(  Stop  ) 


Figure  1.  (Continued) 
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Figure  2.  (Continued) 
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defined  in  Section  I. 

The  three  signals  evaluated  in  this  section  are  a  sim¬ 
ple  2nd  order  AR  process,  and  autoregressive-moving  average 
process,  and  a  sum  of  two  cosines  plus  Gaussian  noise. 

Example  One :  Second  Order  AR  Process 

Example  one  is  a  simple  AR  input  sequence 

x(n)  =  0.75x(n-l)-0.5x(n-2)+n(n)  (64) 

where  n(n)  is  a  Gaussian  distributed  sequence  with  zero 
mean  and  variance  equal  to  one.  This  sequence  should  have  a 
spectral  density  that  resembles  a  second  order  low  pass  sys¬ 
tem. 

A  128  point  sample  of  the  input  sequence  is  provided  in 
Figure  3.  The  expected  results  are  verified  by  the  results 
of  the  Bartlett  estimate,  Figure  4.  Data  from  the  MEM  model 
order  calculations  are  provided  in  Figures  5,6,  and  7.  As 
expected,  AIC  and  FPE  model  estimation  techniques  indicate  a 
minimum  at  order  equal  to  two.  Additionally,  the  Kolmogorov- 
Smirnov  analysis  indicates  a  maximum  at  model  order  equal  to 
two.  These  three  model  order  determining  techniques  clearly 
suggest  that  the  MEM  technique  with  AR  order  equal  to  two  is 
the  most  correct  model.  Figures  8  through  11  show  the  esti¬ 
mated  mean  and  variance  of  the  MEM  and  Papoulis  techniques. 
Notice  how  closely  the  estimated  mean  from  the  MEM  technique 
approximates  the  expected  results  from  the  Bartlett  technique 
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0.00  20.00  40.00  60.00  80.00  100.00  120.00  140.00 

_ TIME  IN  SECONDS _ 

Figure  3.  128  Samples  of  a  Second  Order  AR  Sequence 
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estimation  of  system 


Amplitude  of  FPE  as  a  Function  of  Model  Order  for  Second  Order  AR  Model 


ESTIMATION  OF  SYSTEM 


Figure  6.  Amplitude  of  AIC  as  a  Function  of  Model  Order  for  Second  Order  AR  Model 


Figure  10.  Estimated  Variance  of  Estimated  Spectrum  Using  the 

Technique,  Second  Order  AR  Model 


(approximately  the  true  spectrum).  With  only  64  realizations, 
the  mean  from  the  MEM  technique  appears  to  be  reasonably  accurate 
with  a  variance  depending  on  the  amplitude  of  the  true  spectrum. 
For  this  particular  problem,  the  estimation  errors  are  relative¬ 
ly  small. 

The  Papoulis  technique  fails  to  provide  satisfactory  results, 
with  a  variance  much  greater  than  that  of  the  MEM  technique. 

This  problem  is  attributed  to  the  limitation  of  the  number  of  al¬ 
lowed  iterations  in  the  Papoulis  estimation  routine.  It  is  as¬ 
sumed  that  the  Papoulis  technique  will  do  a  more  creditable  esti¬ 
mation  of  the  spectrum  if  the  routine  is  allowed  to  iterate  be¬ 
yond  100.  This  conclusion  is  suggested  by  Eq(63)  and  Papoulis' 
proof  of  convergence.  However,  since  the  input  in  this  problem 
is  not  deterministic,  the  error  will  not  go  to  zero  as  the  num¬ 
ber  of  iterations  is  increased,  but  will  converge  to  some  con¬ 
stant  value  related  to  the  variance  of  the  noise  process  n(n)  . 

Results  from  execution  of  the  routine  in  Figure  2  are  even 
more  illustrative  of  the  accuracy  of  the  MEM  technique.  Fig¬ 
ures  12  through  15  show  the  results  of  these  techniques  applied 
to  two  input  sequences  of  128  and  32  samples.  The  MEM  technique 
is  clearly  superior  for  this  particular  second  order  AR  sequence. 
Results  of  this  application  indicate  that  the  MEM  technique  es¬ 
timates  the  spectral  density  with  a  high  degree  of  reliability. 
This  can  be  observed  when  the  MEM  analysis  results.  Figures  12 
and  14,  are  compared  to  the  true  spectrum  (by  Bartlett  technique) 
Figure  4. 
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Using  Papoulis  Technique 


ESTIMATED  srEcmun 


Order  Equal  to 


Papoulis  Technique 


Example  Two:  ARMA  Process 

Example  two  is  a  departure  from  the  previous  example,  in 
that  the  MEM  technique  must  approximate  the  system  spectral 
density  [Auto-regressive-Moving  Average  (ARMA)]  by  an  AR  Model. 
This  example  provides  a  more  rigorous  test  of  the  MEM  techni¬ 
que.  Something  other  than  an  AR  model  is  modeled  by  the  MEM 
technique.  Illustrations  of  the  primary  results  are  provided 
in  this  section  with  supplemental  data,  i.e.  AR  model  order, 

MEM  and  Papoulis  mean,  and  MEM  and  Papoulis  variance  estimates, 
provided  in  Appendix  C. 

The  ARMA  sequence  is  generated  by 


x(n)  =  0.5x(n-l)-0.25x(n-2)+0. 125x(n-3)+ 
n(n)-l.ln(n-l)+0.24Ti(n-2) 


(65) 


where  n(n)  is  defined  as  in  example  one.  This  particular 
model  is  the  result  of  a  system  with  three  poles  and  two 
zeros.  It  is  expected  that  the  spectral  data  will  exhibit 
the  results  of  cancellation  of  poles  by  zeros  at  some  frequen¬ 
cy.  A  plot  of  128  samples  of  the  ARMA  process  is  shown  in 
Figure  16.  The  signal  appears  to  be  purely  random.  Applica¬ 
tion  of  the  Bartlett  smoothed  periodogram  indicates  the  ex¬ 
pected  spectral  estimation  results  (Figure  17). 

Application  of  the  MEM  model  order  estimation  techniques 
indicate  an  AR  model  order  of  five  (Appendix  C).  The  MEM 
technique  shows  a  very  close  relationship  between  the  Bartlett 
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128  SAMPLES  OP 


Figure  16.  128  Samples  of  an  ARMA  Sequence  where  the  input  is 

Gaussian  Noise 


ESTIMATE!)  SPECTRUM 


Figure  17.  Bartlett  Smoothed  Periodogram  Spectral  Estimate  of  ARMA  Input 


<>  estimation  and  the  UEM  estimated  mean.  Additionally,  the  MEM 

estimated  variance  is  small  when  compared  to  the  Papoulis  tech¬ 
nique  counterpart.  The  problems  previously  encountered  in  the 
Papoulis  technique  continue  to  be  evident  in  this  example. 

Figures  18  through  21  show  results  of  application  of  a 
fifth  order  AR  MEM  model  and  the  Papoulis  technique  on  128 
and  32  samples  of  the  ARMA  process.  The  MEM  technique  clearly 
estimates  the  spectral  density  with  more  accuracy  than  the 
Papoulis  technique.  The  limitation  of  100  iterations  restricts 
accurate  estimation  by  the  Papoulis  routine. 

Example  Three:  Sum  of  Sinusoids 

Many  spectral  estimation  application  problems  are  a  se¬ 
quence  resulting  from  some  form  of  a  sinusoidal  process.  One 
case  may  be  a  sequence  consisting  of  radar  returns.  This  ex¬ 
ample  exercises  the  spectral  estimation  equations  to  analyze 
a  sinusoidal  time  sequence,  a  process  that  is  the  result  of 
a  sum  of  two  cosines  in  the  presence  of  a  Gaussian  noise. 

This  example  is  set  up  to  provide  an  indication  of  how 
well  the  estimation  techniques  are  able  to  resolve  two  co¬ 
sines  with  different  amplitudes  and  frequencies  that  are 
close  to  each  other.  The  input  sequence  is  generated  by 

x(n)  =  3cos(0.11n)+2.5cos(0.14n)+n(n)  (65) 

where  n(n)  is  Gaussian  distributed,  zero  mean  and  variance 
equal  to  one. 
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ESTIMATED  SPECTRUM 


Figure  18.  Estimated  Spectrum  from  128  Samples  of  ARMA  Process 
Using  MEM  with  AR  Model  Order  Equal  to  Five 


ESTIMATED  SPECTRUM 
PRPOULIS  BANOL1MITED 


Figure  21.  Estimated  Spectrum  from  32  Samples  of  ARMA  Process 

Using  Papoulis  Techniques 


A  plot  of  128  samples  of  the  input  sequence  is  provided 
in  Figure  22.  The  expected  spectral  density  should  have  two 
narrow  spikes  at  0.11  and  0.14  Hertz  riding  on  a  low  level  of 
noise.  The  spike  at  0.11  Hertz  should  have  a  higher  ampli¬ 
tude  than  the  0.14  Hertz  spike.  These  expected  results  are 
verified  by  the  estimated  spectral  plot  from  the  Bartlett 
smoothed  periodogram,  shown  in  Figure  23.  It  appears  that 
the  periodogram  closely  matches  the  spectral  density. 

Estimation  of  the  AR  model  order  indicates  that  the 
best  (minimum  FPE  and  AIC)  order  is  23  (Appendix  C).  Also, 
reasonable  values  for  the  estimated  mean  and  variance  indi¬ 
cate  that  the  MEM  technique  appropriately  models  the  input 
sequence  (Appendix  C). 

The  Papoulis  analysis  indicates  that  this  technique  pro¬ 
vides  increased  accuracy  over  applications  in  the  previous  ex¬ 
amples.  Estimated  mean  and  variance  data  suggest  that  the 
Papoulis  spectral  estimation  technique  produces  reliable  es¬ 
timates  (Appendix  C). 

Results  of  spectral  estimation  of  two  independent  se¬ 
quences,  128  and  32  samples,  are  shown  in  Figures  24  through 
27.  Both  MEM,  with  AR  model  order  23,  and  Papoulis  techniques 
accurately  estimate  the  input  spectral  content.  The  32  sam¬ 
ple  sequence,  when  applied  to  the  FPE  and  AIC  equations,  sug¬ 
gest  an  AR  model  order  of  17.  For  the  32  sample  case,  notice¬ 
able  degradation  in  amplitude  resolution  of  the  MEM  technique 
is  observed.  This  degradation  is  attributed  to  the  large 
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ESTIMATED  SPECTRUM 


Figure  24.  Estimated  Spectrum  from  128  Samples  of  Sum  of  Sinusoids  Process 
Using  MEM  with  AR  Model  Order  Equal  to  23 


ESTIMATED  SPECTRUM 


Process  Using  Papoulis  Techniques 


ESTIMATED  SPECTRUM 


Figure  26.  Estimated  Spectrum  from  32  Samples  of  Sum  of  Sinusoids  Process 
Using  MEM  with  AR  Model  Order  Equal  to  17 


estimated  variance  of  the  MEM  technique  (Appendix  C).  Recall, 
from  Table  I,  the  variance  of  the  estimate  is  inversly  pro¬ 
portional  to  the  number  of  samples.  Therefore,  with  a  larger 
variance  estimate  and  fewer  samples  the  estimated  spectrum  will 
contain  larger  errors. 


Conclusions 


The  signals  analyzed  in  this  section  are  by  no  means  con¬ 
sidered  to  be  an  exhaustive  list  of  all  possible  signals  of 
interest.  However,  enough  data  is  available,  from  analysis  of 
the  provided  examples,  to  draw  realistic  conclusions  about  the 
spectral  estimation  techniques  under  investigation.  The  three 
examples  indicate  the  relative  qualities  of  the  MEM  and  Papou- 
lis  spectral  estimation  techniques.  In  particular,  the  MEM 
techniques  show  considerable  applicability  for  all  types  of 
input  processes  including  simple  AR,  ARMA,  and  sinusoidal  data. 
Tfais  conclusion  is  tempered  by  the  restrictions  placed  on  the 
number  of  iterations  for  the  Papoulis  estimation  technique. 
Clearly,  the  Papoulis  technique  provides  creditable  results 
when  the  input  is  sinusoidal  and  it  is  expected  that  the  tech¬ 
nique  can  be  employed  on  other  types  of  data  if  the  restric¬ 
tions^  are  removed. 

The  analysis  of  short  record  input  processes  is  expanded 
to  specific  problem  signals  from  data  provided  by  the  Rome  Air 
Development  Center.  The  following  section  provides  analysis 
results  showing  the  relative  merits  of  MEM  and  Papoulis  spec¬ 
tral  estimation  techniques. 
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IV.  Analysis  of  Radar  Problems 


The  purpose  of  this  research  is  to  develop  spectral  den¬ 
sity  estimation  algorithms  for  specific  problems  of  interest 
to  the  government.  The  major  limiting  feature  of  the  prob¬ 
lem  is  the  restriction  to  short  sequences  (128  samples  or  less). 
The  specific  problems  evaluated  in  this  section  are  adaptations 
of  data  provided  at  a  recent  RADC  spectral  estimation  workshop 
(Ref  35).  These  problems  represent  data  typically  available 
from  a  radar  problem.  One  problem  is  a  sum  of  sinusoids  with 
a  wide  range  of  amplitudes.  The  other  problem  is  a  sum  of 
sinusoids  with  additive  noise.  Two  different  noise  sources 
are  added  to  the  second  problem.  In  the  first  case,  only 
Gaussian  noise  is  added.  The  problem  is  evaluated  again  with 
additive  Gaussian  noise  plus  impulsive  interference  distribut¬ 
ed  randomly  throughout  the  data. 

Spectral  estimation  techniques  are  required  in  the  radar 
problems  to  determine  the  relative  amplitudes  and  frequencies 
of  the  time  sampled  sinusoidal  data.  The  data  is  assumed  to 
be  defined  over  the  range  of  0.0  to  200.0  Hertz.  One  strong 
sinusoidal  waveform  is  added  to  other  weaker  sinusoids.  The 
total  dynamic  range  is  approximately  40  db  (the  weakest  sig¬ 
nal  is  no  more  than  40  db  down  from  the  strongest  signal). 

Problem  One 

A  simple  sum  of  sinusoids  with  different  amplitudes  is 
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>  '  evaluated  in  this  problem.  The  input  data  is  evaluated  using 

the  interactive  routine  outlined  in  Figure  2.  Independent  se¬ 
quences  of  128  and  32  samples  are  evaluated.  The  sampling 
interval  (At)  is  0.0025  seconds  with  a  maximum  frequency  of 
interest  of  200  Hertz. 

Problem  one  is  generated  by  implementing  the  following  ■ 
equation 

x(n)  =  0.03  cos(2ir78nAt ) 

+0.90  cos(2irl34nAt ) 

+1.20  cos(2irl28nAt ) 

+0.90  cos(2iTl43nAt ) 

+0.04  cos(2irl53nAt ) 

+0.16  cos(2irl65nAt )  ,  n=l,2,...  (66) 

The  dynamic  range  of  amplitudes  for  this  problem  is  -1.43db- 
(-33.47db)=32 . 04db. 

Figure  28  shows  a  128  point  sample  of  the  input  sequence 
x(n).  Numerical  results  from  128  and  32  point  samples  of  the 
input  sequence  are  provided  in  Figures  29  through  32.  The 
higher  amplitude  components  are  clearly  evident  from  both  MEM 
and  Papoulis  spectral  analysis.  Close  examination  of  the 
actual  spectral  density  estimation  data  from  the  MEM  algorithm, 
provide  by  a  computer  listing,  indicate  that  the  MEM  techni¬ 
que  is  capable  of  detecting  the  associated  input  frequencies 
with  few  errors.  In  particular,  the  128  point  sample  clearly 
detects  all  input  frequencies.  However,  MEM  is  deficient  in 


74 


Figure  30.  Spectral  Estimation  from  128  Samples  of  Six 
Sinusoids  Using  Papoulis  Techniques 


ESTIMATED  SPECTRUM 
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Figure  31.  Spectral  Estimation  from  32  Samples  of  Six 
Sinusoids  Using  MEM  Techniques 


PflPOULIS  BflNOLltllTEO 
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Figure  32.  Spectral  Estimation  from  32  Samples  of  Six 
Sinusoids  Using  Papoulis  Techniques 


low  amplitude  estimation,  when  only  32  samples  of  the  input 
sequence  are  available.  The  latter  problem  is  assumed  to  be 
a  result  of  limited  model  order,  i.e.  only  32  samples  of  in¬ 
put,  limits  model  order  to  less  than  32.  A  higher  model  or¬ 
der  may  be  required. 

Papoulis  extrapolation  techniques  estimate  the  spectral 
frequencies  of  the  top  four  amplitude  signals  but  fail  to  de¬ 
tect  the  remaining  signals.  The  Papoulis  technique  becomes 
less  reliable  with  fewer  input  samples,  a  problem  encountered 
in  the  last  section  and  related  to  limitations  imposed  on  the 
number  of  iterations. 

An  additional  note  of  interest  is  the  failure  of  MEM  es¬ 
timation  to  estimate  the  relative  amplitudes  of  the  various 
frequencies.  This  problem  is  clearly  evident  in  this  and  pre¬ 
vious  examples.  Ulrych  and  Bishop  (Ref  39:191-192)  attribute 
this  anomaly  to  the  relatively  high  variance  of  the  MEM  spec¬ 
tral  estimation  techniques  when  compared  to  other  techniques 
(also  see  Table  I).  The  variance  is  directly  proportional  to 
the  amplitude  of  the  estimate.  Ulrych  and  Bishop  suggest  that 
this  problem  is  allowable  if  a  trade  off  between  variance  and 
resolution  is  considered.  Clearly,  the  results  suggested  in 
Section  II  are  illustrated  in  this  sum  of  six  sinusoids  prob¬ 
lem.  If  the  Papoulis  technique  successfully  estimates  the  in¬ 
put  frequencies,  the  relative  amplitudes  of  the  generating 
data  is  preserved.  Analysis  results  from  Section  III  and  this 
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section  indicate  resolution  limitations  of  the  Papoulis  tech¬ 
nique  when  compared  to  the  MEM  technique. 

Problem  Two 

The  second  radar  problem  considers  the  case  of  two  close¬ 
ly  spaced  cosines  with  additive  noise.  The  second  part  of 
this  problem  includes  an  additional  random  noise  source  that 
is  impulsive.  Initially,  the  sequence  of  interest  is  generat¬ 
ed  by 

x(n)  =  231  cos(2ir51nAt ) 

+1314  cos(2ir62nAt)+n(n)  (67) 

where  n(n)  is  a  Gaussian  process  with  zero  mean  and  variance 
equal  to  36.  The  large  variance  indicates  a  large  noise  source. 
Figure  33  shows  a  128  point  sample  of  the  input  sequence.  Spec¬ 
tral  estimation  results  from  128  and  32  point  input  sequences 
are  provided  in  Figures  34  through  37.  Satisfactory  results 
are  indicated  for  both  the  MEM  and  Papoulis  techniques  applied 
to  the  128  point  sample.  The  Papoulis  technique  fails  to  dif¬ 
ferentiate  between  the  two  input  signals  when  only  32  points 
are  analyzed.  This  problem  is  another  indication  of  the  resol¬ 
ution  differences  between  MEM  and  Papoulis  spectral  estimation 
techniques. 

An  additional  noise  process  is  added  to  Eq(67).  Samples 
from  a  uniform  random  number  generator  are  taken  at  each  point 
n  .  The  uniform  generator  produces  numbers  between  [0,l]  . 
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Figure  35.  Spectral  Estimation  from  128  Samples  of  Two  Sinusoids 
Plus  Noise  Using  Papoulis  Techniques 
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Figure  36.  Spectral  Estimation  from  32  Samples  of  Two  Sinasoids 

Plus  Noise  Using  MEM  Techniques 


Whenever  a  number  is  greater  than  an  arbitrary  value  of  0.85, 
an  additional  value  of  9999  is  added  to  the  output  x(n) 

This  impulsive  signal  is  added  to  simulate  external  pertur¬ 
bations  such  as  lightning  and  erroneous  data  coding.  It  turns 
out,  that  an  average  of  ten  pulses  occur  in  every  128  samples. 
A  128  point  sample  of  this  sequence  is  shown  in  Figure  38. 

The  data  are  evaluated  using  MEM  and  Papoulis  spectral  esti¬ 
mation  techniques.  Figures  39  and  40  provide  graphic  examples 
of  the  resulting  analysis.  Both  techniques  are  unable  to  re¬ 
solve  both  frequencies.  However,  MEM  does  indicate  a  large 
amplitude  node  near  the  two  frequencies  of  interest,  i.e.  51 
and  62  Hertz.  Note  the  very  large  impulse  near  zero.  This 
added  impulsive  noise  reduces  the  signal  to  noise  ratio  of 
both  techniques.  Therefore,  the  resolution  of  MEM  and  Pap¬ 
oulis  techniques  is  greatly  reduced. 

Conclusions 

The  problems  considered  in  this  section  represent  "real 
world"  radar  applications.  Papoulis  Bandlimited  Extrapolation 
and  MEM  are  capable  of  resolving  many  of  the  signals  present 
in  the  input  sequences.  In  particular,  a  sequence  of  128  sam¬ 
ples  can  be  resolved  when  the  maximum  deviation  between  fre¬ 
quency  amplitudes  is  less  than  40  db.  Additionally,  when  res¬ 
olution  is  of  greater  importance  than  amplitude,  the  MEM  tech¬ 
nique  exhibits  marked  superiority  over  the  Papoulis  algorithm. 

Noise  perturb  'tions  do  degrade  the  relative  accuracy  of 
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Figure  38.  128  Samples  of  Sum  of  Two  Sinusoids  Plus  Noise 

Plus  Impulsive  Noise 


ESTIMATE!)  SPECTRUM 
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Spectral  Estimation  from  128  Samples  of  Two  Sinusoids 
s  Noise  Plus  Impulsive  Noise  Using  MEM  Techniques 


EST MATEO  SPECTRUM 


Figure  40.  Spectral  Estimation  from  128  Samples  of  Two  Sinusoids 
Plus  Noise  Plus  Impulsive  Noise  Using  Papoulis  Techniques 


each  estimation  technique.  In  particular,  considerable 
resolution  discrepancies  are  observed  when  both  impulsive 
and  normal  noise  are  introduced  as  compared  to  simple  normal¬ 
ly  distributed  noise. 
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V.  Conclusions  and  Recommendations 

Power  spectral  estimation  of  short  record  time  series 
data  represents  a  unique  problem.  Traditional  estimation 
techniques  such  as  FFT  and  periodograms  do  not  provide  suf¬ 
ficient  resolution  for  adequate  analysis.  Two  techniques  are 
suggested  that  produce  spectral  density  estimates  that  are 
clearly  advantageous  and  computationally  desirable.  The 
problems  considered  in  this  research  are  limited  to  station¬ 
ary  processes.  The  Maximum  Entropy  Method  applied  by  Burg 
and  bandlimited  extrapolation  by  Papoulis  clearly  are  desir¬ 
able  algorithms.  Each  technique  exhibit  some  limitations 
and  certain  trade  offs  between  resolution  and  amplitude. 

The  spectral  estimation  techniques  are  applied  to  sev¬ 
eral  "real  world"  problems  producing  acceptable  results  in 
most  cases.  The  Papoulis  techniques  do  not  provide  the 
resolution  of  the  MEM  technique,  but  this  reduced  resolu¬ 
tion  is  attributed  to  the  prelimiting  of  computer  process¬ 
ing  time. 

This  research  does  not  suggest  that  all  possible  analy¬ 
sis  techniques  are  considered.  Several  techniques  are 
currently  being  evaluated  by  several  individuals.  Of  par¬ 
ticular  importance  is  the  work  by  Makhoul  et  al  (Ref  26)  in 
the  area  of  lattice  networks. 

-«  Additional  areas  of  evaluation  and  research  are  extensions 
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to  methods  covered  in  this  thesis.  Extension  of  the  techni¬ 


ques  to  complex  signals  should  be  a  relatively  simple  operation 
Also,  since  signals  from  imaging  systems  can  be  evaluated  using 
the  techniques  of  this  thesis,  two  dimensional  spectral  esti¬ 
mation  techniques  are  a  reasonable  evolution.  Finally,  Cadzow 
(Ref  14)  suggests  a  technique  to  evaluate  an  input  sequence 
using  MEM  when  a  percentage  of  the  data  is  incorrectly  coded. 
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Appendix  A 


Prolate  Spheroidal  Expansion 


Prolate  Spheroidal  expansion  methods  can  be  used  to  de¬ 
termine  P(u>)  from  a  sample  of  an  input  sequence  xx(t) 
xi(t)  can  be  expanded  into  a  series  of  values  dependent  up¬ 
on  a  set  of  orthogonal  eigenfunctions  <j>k(t) 


a,  = 


1_ 

X, 


xi(t)$k(t)dt 


(68) 


where  are  the  eigenvalues  and  T  is  the  sampling  inter¬ 

val  defined  for  the  Papoulis  extrapolation  in  Section  II. 
xj(t)  is  expanded  as 

00 

Xl<t)  =  XVk(t)  (69) 

k=0 

since  4»k(t)  are  orthogonal  eigenfunctions,  the  Fourier 
Transforms  are  a  scaled  value  of  <J>k 

♦k(t)*  4>k(bu)Prt(a))  (70) 

k  /T~  k  a 

/Ak 

B  =  /Stt'T/o  ,b*T/a 
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(71) 
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is  the  energy  in  x(t)  and  is  finite.  Since  the  eigenvalues 
are  less  than  one,  1-Xk<l  and  1-Xk^l-XN  when  k£N 
Therefore,  the  energy  in  the  nth  iteration  of  the  Papoulis 
technique  is 

En  *  I 

k=0  k>N 

N 

<(1-XN,!“  I  “k"  *  Z  “kZ<e  <76> 

k=0  k>N 

where  e  is  some  arbitrary  value.  Therefore  (l-XN)2n-*-0  as 
n-M<>  and  Pn(u>)  converges  to  P(oj)  as  the  number  of  itera¬ 
tions  increases. 
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Appendix  B 

Subroutines  Used  for  Spectral  Estimation 

Several  algorithms  are  used  in  the  analysis  sections  of 
this  paper.  This  appendix  provides  a  brief  explanation  of 
each  routine  along  with  a  FORTRAN  IV  listing.  All  subroutines 
are  compiled  into  object  code  and  stored  in  an  ASD  computer 
system  program  library  called  MEMBIN.  A  listing  of  the  sub¬ 
routine  names  and  functions  is  provided  in  Table  II.  The  ap¬ 
plicable  equations  for  each  routine  are  repeated  in  this  sec¬ 
tion  from  section  II. 

Subroutine  AAIC 

Calculates  the  Akaike  A-Information  Criterion  using  Eq(54) 


(AIC)p  “  N  In  (Pp+1)+2(1+P)  (54) 

where  N  is  the  number  of  samples; _ln  is  the  natural  logrithm, 
p  is  the  particular  model  order  and  Pp+1  is  the  innovations 
(estimation  error)  power. 

C 

SUBROUTINE.  AAIC(N,h,P,A2  C) 

c 

C  AmIC  CALCUwfcTFS  THE  AKAI<t  A  •lNF0RH*t  fl  lfl  CRITERION 
C  VALUE  FROM  INNOVATION  EfrCK  VARIANCE  (POWER) 

C  VECTOR  P()  (1,2). 

C 

C  CALL  c.0  BY  1  CALL  AAIC  (  N,  H , P, A1C) 

c 

C  N  IS  THE  1 D T 6 L  NUrtBEF  UF  lHf-UI  SAMPLES 
C  M  IS  THE  •ODcL  URGE®  PLUS  Oi- E 

C  PI)  ARE  THr  FrPCR  P-JWCK  VALUS  FOR  OK  )E"*C  U*>  TC  M 
C  AICO  AI..E  r  HE  CALCULATED  INFORMATION  CR3  T  CRT  OK  VALUES 
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TABLE  II 


Subroutine  Identification 


Ex 

Routine  Name  Function  Performed 

ternal  Calls 
Required 

AAIC 

Calculates  the  Akaike  A-Inform- 
ation  Criteria  (ATC) 

None 

AFPE 

Calculates  the  Akaike  Final 
Predictor  Error  (FPE) 

None 

BTLET 

Computes  a  Bartlett  Smoothed 
Periodogram 

FFT 

FFT 

Decimation  in  Frequency  Dis¬ 
crete  Fourier  Transform 

None 

HGRAPH 

Generates  plotting  Data  for 
CALCOMP  x,y  data 

AXIS,  SYMBOL 
PLOT,  SCALE 

KOLSMR 

Calculates  the  Probability 
that  the  MEM  Coefficients 
Produce  a  True  estimate. 

PDF , NKS 1 
(IMSL  Routine) 

MEMC01 

Calculates  the  Innovations 

Power,  Estimation  Coefficients, 
and  Autocorrelation  values  for 
the  MEM  Technique 

None 

PAPO 

Calculates  the  Spectrum  using 
the  Papoulis  Bandlimited  Ex¬ 
trapolation  Routine 

FFT 

PDF 

Calculates  the  Values  for  a 
Theoretical  Uniform  Distribu¬ 
tion. 

None 

PSPEC 

Calculates  the  Power  Spectral 
Density  using  the  MEM  Estima¬ 
tion  Coefficients 

None 

10 


2 


*ii  j  riM  15  BIST  9BAI.ITT 

t&M  0UTT  fi55is»*3  TO  DDC  . 


C 

C  1.  4<AI<i,J.fn4  NEW  LOOK  kT  THE  S^Ti^Ti CAL  MODEL  IDEN7IFI- 
C  AT  a  0N"»  ItTE  TRANS.  Al’TCM.  COhT.»VJL.<  C-l  ■},  r  90-7 »  PEC  *’  ‘  . 

C  2.  JuNEs,  ?.w.»  "lOEtr  T“iC»'-Ti  CN  AN  3  A  J*i  Of  r“’cSrI  VE  SPtCT<UM 
C  EST  I  :A1  rf"»",iEL:Z  l.-vANS.  #-UTOf.  CO  U  •  ,  VdL  •  AC-19.  ?  if, -23  . 

C  DEC 
C 

c 

DIMENSION  P(l)  ,P  IC(i) 

C  G E WE »\ T  E  T-l"  ^TANTiUG  i(«ut  F„r<  iH 
AlC(i)sN*'ALOG(o(i)) 

IF  (H.E3.  i  )  KET'J'N 
C  GET  THE  4jC  VALUES 
DO  l j  i=?,M 

AiC(i>=S‘ALUG(P(f  ))  +2->  (i  -1) 

13  CD.iriJJF 
t-cTJ  .\N 
END 


.Subroutine  AFPE 


Calculates  the  Akaike  Final  Predictor  Error  using  Eq(52) 


(FPE), 


Btt  Vi  •  0  P*N/2 


N+P+l  (-  Vi  ^  N/2  pin 

[N-P-l  \^(1.5-P/N)y  »  N/Z  PSN 


(52) 


where  N  is  the  number  of  samples,  p  is  the  particular 
model  order,  and  ^p+i  is  the  innovations  power.  This  e- 
quation  is  an  adjustment  to  Akaike's  FPE  formula,  suggest- 
ed  by  Jones. 

Q  +  + »*+*,+'-*••*'+*•++*++  +  +*  *******  *■*  +  ***  ■*>*<■***  ♦  *♦«♦»♦♦»«** 

C 

SUBROUTINE  AF?£(r  (i,P,FFE) 

C 

C  AFP£  CALCULATES  THE  AKAIKE  FINAl  PREDICTOR  EFFOR 
C  CALLED  ?Y!  PALL  AFCfc  (N,;1 ,  P,  FFE) 

C 

C  N  IS  THE  T  D T  A  L  NUMBER  CF  SAFFLE.» 


Mil  fin  IS  BIST  All  ITT  PTUilKA 
man  earT  mbi»is>*«;;  tv  wc 

c  rt  is  the  »:d^ci.  order  plus  o?  i 

C  P()  ARE  T,JE  "PkDR  T„\  r>’/ATICf  ?  (POM^K) 

C  FPEO  ARE  T  MP  CALCULI!  tC  FPI  VALUE1? 

c 

C  1.  AKA  IKE*  4  •  » “FITTING  oUTCrT  GiEobiVE  MODELS  FT  R  PF  E0lCTl  OfT » 
C  ANN.  INS'*.  S 7 A T  .  - at  w.  ,  VC  L  . 2i  ♦  2  7-7  ,3  -V  ?. 

C  2.  JouS<,  ?.«.*“AUT&*  cGi-.^SSiVL  3ZL  EOT  10!  ",  GE0PHY3.  , 

C  VLL.-.i, :  /1-7.AU3  •  •  . 

c 

c 

01  •*Erl37f'\  P(l)  »rc  £ ( 1 ) 

C  GENERATE  T-l?  ^'CALE  F6CTCR 

FPE(l)r(FLyAT  (N  + 1  )/Fl«j*.T  <K-i))*P(  1) 

IF  RETURN 

C  NORMALISE  *rFULTS 
FtE’<p=-pf  (1) 

FP£  ( 1)  =1  • 

C  GET  T^r  rs-  VALUES 
NNsN/2+i 
OC  i:.  I  *  2  f  M 

FP£(I)= (FLOaT (N+I-D/FlCAT (N-I-l) ) *P( i) 

IM 1=1-1 

IF  (iHi.  L*  «NN)  GO  T0  Z2 
FPt(I)=rPE(I)/(l.;;-(I-l)  /N) 

20  FPE( I ) =F°r ( D /FT  t«D 

10  CONTINUE 

N£  TURN 
END 

Subroutine  BTLET 

BTLET  determines  the  estimated  spectrum  by  applying  the 
Bartlett  smoothed  Periodogram  equations. 

V“>  ■  s  iv->r 

K 

Vw)  b  S  E  ^pi(w) 

i«=l 


where  X^w)  is  the  N-point  DFT.  There  are  k,  N-point  per- 

A> 

iodograms,  Ppi^) 


(11) 

(12) 
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.***■»•  *.  ***♦*.  .«♦.*♦.  .i*.-  *■»♦.**■*  *•-•*»<. 


SUVUJTINE  RTLET  <M,F,  i, or  N, 31  ,flo,fm:  ,x,  » 

c 

C  37lET  C’lM^iir  e<?  Thf  ES’IFMtD  FOWE''  Sn?TR'.H  USING  THE 
C  BART LET  t  SMOOTHED  °l!  I  DT'OCrA  P  *ETHOO  <11. 

c 

C  CALLED  BY*  C.'LL  3TL  E  1  < N,  .•:,»  C  CM,  OT,  FLO,  CHI ,  X, N  ) 

C 

C  N  IS  1 H-  L"»MPER  3  F  SAMPLES  PE*  WINDOW 
CM  IS  T'U  ‘"IMBfc-  OF  PF  J.  I  Of  C  GR#;  MS  USED 
C  10W  IS  IISFO  ~Q  COMPUTE  N  Fin  RJUUNE  FF  T  ,  Ns?**HCDN 
C  DT  13  TM;  ttme  SAHCLc  IuTcKVAL 
C  FLO  IS  THF  low  Fv-our\fY  OF  a.^7  E *E6f 
C  FHI  13  TH?  HIGH  FaE  CUE  N'C  Y  IF  iNTERE  5T 
C  XO  ARC  T“-  f*FQ'JFNCV  SIKH  FS 
C  Y<)  ATE  THE  fiHFLI  PIPES  A  7  *»H£  FREQUENCE  E> 

c 

C  KcOUInEO  i  O'JTTKjtS  P.-C^  Llr-r.hFYt  FF7 
C  REQUI  .  ES  <\  +  -  )M  DATA  PUNTS  IM  LOCAL  F:  L“  ”  A  FEb  • 

C 

C  1.  OPPEUH-:r^,  A.v.  A  *  D  W.  SC  H^FEa,'*  DIGITAL  SIGNAL  PPOIES 
C  ING*V4cU  VOnKt  Pnc  EM  1C  c-F  ALL  t  13jL  . 

c 

*<*-  ■» fc j  *.».  ■» -  *.  ,  ».*v  *■  •••  »•  ****»■  ♦ 

C 

DIMENSION  X  (1)  »  Y  ( 1) 

COMPLEX  P,XF<3*.) 

tiP  j  j  =N  ♦  I  0 

MMPl=MOOM+l 

NM*2*tf 

NP1=N*1 

C  OF  IS  THE  CHANGE  I '4  Fr.EOUENfY  FOR  EACH  S AM°LE. 

OF* (FHI -FLO) /FLOAT <N> 

FREQ*. . 

oo  n  i*i, n 

Yd)*:. 

X (I) =r REO 

10  FREO=F\EO+OF 

00  22  I* 1 , M 

00  3*5  j  =  i,np»:. 

READ  <c  ,  3  0  )  XI  N 
IF  (J«Lt.  Of  )  GO  TO  3C 
XF(J-.jf  )  *XIN 
30  CONTINUE 

00  1>  J*M°1 ,NS 
15  XF(J)=C.,C.) 

C  GET  THE  OFT 

CALL  FFT <XF,MU°i,i.) 
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fr:) .  «  "Ut  lSSSST  WALlTt  ntiU>fTJ 

00  “  J  =  1  ,  N  ** 1  *°  300  ' - ^ 

«=C0fJJ?(yc-(  j) ) 

A0V=  Xc  {  J)  A 
40  Y  (  J)  =Y  ( J )  +a  3V 
2  c  CONTINUE 

DO  .>c  I  =  1 ,  N 

50  Y(I)=Y  (I)  /  (FLOAT  ID-FLOP  1  (N)  ) 
wLWINQ  s 

600  roi«!AT(P2'  .l«!) 

k:  t*j -;im 

£40  ..  .. 

Subroutine  FFT 

FFT  calculates  a  2M  point  DFT  or  IDFT  of  the  input  data. 
The  routine  is  a  standard  decimation  in  frequency  application 

C ***  **■  *  ■  '  f>  f  '  i  •>'  ^  0  M  .  ^  ^  ^  4,  •  «  4  4  •  I  •'  k  O  |ir.  i  a  4  |  A  Jl  ^  a  4t 

C 

S'J3'OU::‘4£  FF’(X,m,  Xi  J 

c 

C  FFT  C0HDUT=S  "HE  DPT  FOR  a  VECTOR  INPIJ”  X()  T»aT  HAS 
C  LCGOASF  ?)  M  POINTS.  PCr  FOiW/RD  FFT ,  XI  =  i.,  FOP 
C  REVERSE  (TOFT)  XI  =  -1. 

C  THE  ROUTINE  IF  h  DEC  I  HAT  I  Oh  IN  FREQUENCY  (i>. 

C 

C  CALLED  3Y!  CALL  FFT(X,R,XI) 

C 

C  XI)  ARE  THr  INPUT  VALUES  TO  EE  TRANSFORMED,  RESULT 
C  IS  IN  XO  ON  RETURN 

CM  IS  THE  MULTIPLYING  FACTOR,  M=LuS(OfSE  2)4' 

C  XI  DETERMINES  IF  A  FFT  Or  I FFT  IS  PfRFCF  TED, 

C  XI  =  1.  IS  FOR  FFT,  XI  =  -1.  IS  FOR  TO FT 

C 

C  A.  UPPEMHETM,  A.v.  A  NO  H.W.  SCHAFER, "DIGITAL  SIGNAL  PROIFSS 
C  1 NG",  NEW  YOkXJ  PRENTICE-HALL,  i>l?h. 

c 

4  4,1  r  .,.1  4.  f  •  ,  41  +  44*  +  *■  ♦  * 

c 

COMPLEX  X(1),U,H,T 

N=2‘*M 

NV2=N/2 

(nM1=N-1 

C  BIT  REVERSAL  PART 
J=1 

PI*3 .1  '■*  1  r9  2CE  3Kfl  °£» 
no  7  1=1, NM1 
IF  (I  »GE,  J)  GO  TO  l 
T=X  <  J) 

x(j)  =xm 

X(I)=T 
5  K=NV  2 
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fcsa  1 3  I  ■  - 1 £ 


LI  5JSI  nuyiW#*** 

KaXlaHSD  Tv>  •- — 


6 

IF (<«GE« J) GO  TO  ? 

JsJ-K 

KsK/2 

GO  TO  ; 

7 

J=  J  +  K 

CO  2‘.  L  =  1,M 

LES2  T ' L 

Lti=LE/? 

U= (!.,>.) 

WsCEXp(CH3LX(  ,.,-XIJ  PI/L  E.  1 ) ) 

c  atm 

t'  FL Y  CALCULATIONS  FA h  T 

00  Zn  J=  1 » L  El 

QO  1*  t=j,n,le 

I®=I+LEl 

T=X  (I3)  •  1) 

X  ( Ip)  =  X ( T) -1 

1C 

X(I>=*<I) *T 

20 

»J=’J  w 

iF(XI.3T.  ,)RET!)cw 

DC  30  I=i,N 

30 

X ( I ) =X( T) /N 

Me TURN 

END 

Subroutine  HGRAPH 

Given  some  input  data,  x  and  y,  and  the  appropriate  lab- 
ling  information,  HGRAPH  provides  an  8J  inch  by  11  inch  CAL- 
COMP  plot  of  the  data  (Ref  42) 

■»*«.#  »««•««  »<**»*•  I.  v-»*  *  4  v.  <*•*♦*»■**-•  *■  4***-«..*4*»** 

C 

SUa^OUTlMS  HGRAPH(X,Y,M,  iO»NO,NP*N3) 

c 

C  HGRAPH  P-UVinrs  A  Cm.-C0;'°  dLCT  uF  THE  IMP  I T *3  Y  ANO  Y 
C  WITH  7H=  L*RL"S  e.'.OVIQEP  =?Y  *rRAY  ID  (i). 

c 

C  CALLEO  -lYt  CALL  HGRA  PH  (X  ,Y,  N  ,iU  ,HJ.NP,Nr,) 

C 

C  HGRAPH  M'lET  R"  PRECEEPEO  EY  CALL  PLOT  -3)  AND 

C  CALL  »LCT((». 03,-3)  INITIALLY  TO  GET  UP  CALCOMP  PLOTTER. 

C 
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0-1.3  IS  BIST 


XO  a\f  TH“  X  OhOi^er-FS 
Y ( )  ARE  THf  Y  ORDINATES 
N  IS  THE  NUM^'P.  OP  PliNTS  In  X,Y 
IOO  ARE  TM"  l  4 *LE  DATA  V ALL  E  S 
NO, N»,  VS  ARF  PEN  CONTROLS 

REOU I\SO  73'JTTNESI  SCALE,  PLt  1  fSYM30L»LI'ir.* 

1.  VAUGHN,  3.  ,  LECTURE  N  A  T  £  k  I A  L  £E6.  9 1 , 0 IG ITAL  SIGNAL  PRO 
E3SIHG,  cr:MOOL  OF  ENGINEERING,  AFIT,  lA?f>. 

Q4  *  *  »  ♦  -,  »  +  >  4  '  «  '  !  »  *  4-4"*  A  4  ..  4.  4  *  4  4  4  *!-»»»♦*  *••*.,  ~  -*.«.«  .  4444 

Ol'-r .NSION  X  (1 )  ,Y  (1)  ,IL(J  ) 

if (NO, so.?) CALL  PLOT (-1.6,, 2.1 3) 

IF(NO.Sn.?)G0  TO  3? 

IF(:(0.LT. ')  GO  "0  1. 

cm.l  s:ale(x,?.,n,d 

CAuL  EGA  Lr(  Y,0.,N',1) 

10  Cm  LL  =>LOt  <  C  .,  11.  ,  2) 

C'llL  PL0T  (5.S.,11.,Z) 

CALL  PLOT  (a  .,  ,  0.  ,  2) 

CALL  PLOT Cfi.tP.,?) 

CAuL  PL  OT  (1  •  3r’ ,  1  •  3:>  , *  3) 

CALL  PL  OT  (  L  .  ,  6 . 3  L  ,  -?) 

IF  (10(1)  .  EQ.9S9)  Go  TO  2: 

CALL  PL  OT ( . 1 , - . 1 , -3 ) 

CALL  PLOT (0  .,-2. ,-2) 

DO  2T  1  =  1, 7, 2 

20  CALL  SY'-fOOL  <<1+1.  :)  *  .  1 ,  .  3  , .  - 7  , ID  ( 1 )  ,  90  .  ,  2f  > 
call  plot ( o  •  *  c  • ,  2 ) 

CALL  P..0T  (l.,<  .,2) 

CALL  PLOT(l.,2.,?) 

CALL  PwOT(.l.,2.,-2) 

CALL  eLOT(-.i, .1,-3) 

2&  CALL  PLOT  (3.6,-.).  , -2) 

CALL  PLOT(C  3  ,-2) 

CALL  PLOT (-5.8, j. , -2) 

CALL  5YNflOL(#t,-.2,.l,:r  (13)  ,  j.,5D 
CALL  PLOT(&.3,.7V  ,-3) 

CALL  4XrS(3.,(..f  1D(V)  ,-2i  ,7  .  ,«K  .  ,  X  (r>  *1)  ,  Y  (N +  2) ) 

CALL  4XTS(0.,".*:0(11)  ,2ii,i;.,18:.,Y(i4+l)  ,  Y  (N+2)  ) 

30  Y (N>2) =-Y (N+2) 

CALL  LiHF(Y,X,N,l,HF,NS) 

Y (N  +  2)  =-Y  (N*2) 

CALL  °LPT(i.c.>,-2.1t,-3> 

RE TURN 
ENO 


ftvu.tn  matm** 

«cn  TO  DOC  •  — p- 
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11 


Subroutine  KOLSMR 


Given  a  set  of  M  estimation  Coefficients,  determine  if 
the  AR  model  accurately  models  the  input  data  sequence.  An 
error  vector  is  generated 

e(n)  =  x(n)-x(n)  l^nlN  (48) 

where 

x(n)  =  £  ®pJx<n-3>  <47> 

j-1 

Assuming  Eq(17)  is  correct  we  proceed.  If  x(n)  accurately 
estimates  the  input,  e(n)  will  be  a  white  noise  process. 

The  routine  generates  e(n)  then  calculates  the  spectrum  by 
the  Bartlett  smoothed  periodogram  technique.  A  bandlimited 
sample  of  the  frequency  data  is  integrated  (Simpsons  integra¬ 
tion)  to  produce  a  quasi  distribution  function.  The  distri¬ 
bution  function  is  evaluated  by  a  Xol  mogorov-Smirnov  two 
sided  test  to  determine  how  well  it  approximates  a  uniform 
distribution.  The  results  are  a  probability  measure  of  good¬ 
ness  of  fit. 

c 

SUBROUTINE  KOLSMR  ,  F7A  ,P".3  *> 

C 

C  KOLSKs  C 0  1CUTES  THE  cR'Or'A  BiL  ITY  THAT  THE  SYSTM  OF  uRO"? 

C  M  IS  CO^EOT.  THE  TEST  IS  ACCCHPLISHEO  OY  TH*  KUL  MOGOBOV- 
C  SMIRNOV  SINGLE  SA  *iDL  £  1ES1  F  OF  G00QM255  PI’7  (1). 

c 

C  CALLED  BY?  CALL  KOLSH*  (HCDf  ,  M,  ST  A,  PR9M) 

C 

C  (N+5C)5j  POINTS  OF  DATA  MUST  CE  AVAILABLE  TN  A  LOCAL  FILE 
C  TAPE f • 


•WiS  TL&t  IS  BIST  aVU.ITT 

«»ri  FtfAAiittJti  »  jUfl - 
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c 

C  MOON  I?  T-»c  *  Tc  CAlCUlATE  it  OF  SA  H=>L  Er 
CM  13  r+E  ‘•‘ODEL  O'DLr.  FL'iS  ONE 

C  8TAO  ARE  THF  ESTIMATION  COEFFICIENTS 
C  PkBM  IS  3%0rA9I  L  IT  Y  1  HA  1  THE  SYSTEM  MOO'D  IS 

C  CORRECT 

c 

C  ROUTINES  PE^UIREOt  Pr:p,  NK  Si  (  AN  I  MSL  ROUTINE) 

C 

C  1.  30X,  O.E.  A  NO  1.  JENKI?-S,»TInE  SERIES  ANALYSIS  FGkCA?THC- 
C  A  NO  CO.*Ti  OL**»HOLOt  N-CA  Y ,3  S“u» 

C 

(J**4  +  *4  •  «•  '  ♦  *•••»  r.  *i  *  «  -.  ‘  r  4v  il  ■  W4«  *-  A  ►  * 

C 

C  PDF  IS  CALICO  9Y  JHSl  kOUTINc  NKSl,  THE  KOLO-CMlR  kGUTTME. 
EXTERNAL  nO  F 

Oi  HENS  I  i.'M  yv(lJr  )  ,E.'  (  i«,  )  »01  A(l)  *E(i  IjTJ),PD1F(o) 

COMPL-V  A,XF(3C-A) 

N=2*  %HODN 

C  SKIP  .  •]  "OTUTS  OEF'J-  t  EACH  H I N'OO W 

NP->.,  =  ^4r;'. 

C  ZERO  OUT  THE  ARRAYS 
DO  1*1  1  =  1,  M 
10  XM ( I ) s  C  # 

DO  2 A  J  =  1,N 
XF  ( J  )  =  (  '  .,j.) 

20  fcR( J) =" « 

C  CALCULATE  THE  ESTIMATED  OUTFUT 
C  EST  OF  X(I)=SIIM  (PTA  ( I )’  X(l-1)>  (I  =  1,M) 

X=J. 

M<«l=M-i 
MM 2= MM 1-1 
NP3jX>  =NP«|?ii*  Sf- 

C  WE  WANT  3-'  WINDOWS  OF  N  SwMPLLS  1'z.rs  IRS)  E A°H 
DO  ?u  1  =  1, NP*.  :xc 
iFMMl.LT.2)GO  TO  nf 
no  A:;.  J=1,MM2 

40  XM(MM1-J*1) =XM(HM1-J) 

45  XH (1 )  =  X 

SUH= j. 

00  J=t,MMl 

50  E'UM=  3UM  +  PTA  ( J  +  l)  •  XM  (J  ) 

Rf  AD  (c  ,  >  n  • )  X 

C  FIND  THE  ERROF  FO<n  THIS  SAMPLE 
30  E(I)=X-StiM 

C  NOW  CO  TH=  9 A F'T  LETT  LhOOTHEC  FEh  1000SRAH  OM  FT  WINDOWS 
C  N  POiN  rs/wiNnow 
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■;W* 


IS  W**1** 


tuxu*** 


TO 


-I*1’ 


av.j. 

00  *=1,5...  .**,**  »'V' 

ou  J=i,NP;:, 

IF  ( J  » L  5!  #  •>  >G0  TO  To  fP-'> 

xf  ( j  -: '  >  =  :  <  j4  <  * p  ;  i  -m  pi  t ) 

3  CO  NT  I  >-"..•=■ 

CTLL  FT" (yf,M00M,1.) 

DJ  d '!  J=i,N 
A=CONJG( YF(J) ) 

A«3V=XF(J>-  A 
3  £R<  J)="R<  I)  +A9V 

;  continue 

oo  j-  t=i,n 
£P(I)  =  Fr-(T)/(h  ..  N) 
i  AV  =  i  V+E~ ( 1 ) 

WE  NOW  S  SMOOTH  c[!  SPEC.  kUM  Oc  THE  Fr^P,  IF  WGN, 

WILL  PE  rL&T.  LET  Tu-  ,,  y/^o.L’LS  ..tPKLScNT  A  OFT  SITY 
FUNCTION  F O'v  SUM  < i,\J. 

SCALE  AND  I N T  rG  i;  t  T  S  "0  C-tlT  T  HE  IS7:-.I  QU”  TON  Or  THE 
PKJC ESS. 

SUM=*. 


O'-  THE 


CF;Mi.=  .  . 

F\M£  =  . 

C  00  A  SIMPSONS  INTEGRATION 
00  11.  1=1, N 
£KC)8«(I)/AV 

SUM=S'JM+(  ERM2+4.  '  ERMl+E^  (D)  /6. 

ERM2=S\M1 

ERHls'-Pd) 

110  £*.<i)  =  S!IM 

C  NOW  1 EST  F0»  UNIFORM  OVER  THE  INTERVAL  (  *  TO  l.G  ) 

CALL  N<S1 (PDF, £<,N, PDIF,  1LR) 

PRPM-PniFtb) 

REW1NT  5 

*  5Cu  FORMAT*  r  2".  It  ) 

RETURN 

END 


Subroutine  MEMC01 

With  an  N  point  input  sequence  of  data  and  some  assumed 
model  order  P  ,  determine  the  estimation  coefficients  §  ^  , 

innovations  power  Pp+1  »  and  the  autocorrelation  coeffic¬ 
ients  4>A  ,  this  routine  carries  out  the  recursions  outlined 

by  Anderson  (Ref  6). 


Ill 


C *♦  4  »  ♦  "  '  <  f  *  *  -  ♦  M  •*  j.  4  4  ■  *  *. .  •»  -I  -.  - 


*  4  >*■'»■•  » 


f  »  ■»*■  *  *4  V*  ♦  -•*» 


c 

SUBROUTINE  Kr.'iC 01  (N» I1,X,  FHI  ,CTA,P) 

C 

C  MEMO  01  CALCULATES  fHF  £SIIK;T10:i  VALUES  FOR  A"  MTH 
C  ORDS  r’  INPUT,  MCMCCt  IS  Al  ADAPTATION  of  am  algorithm 
C  BY  N.  AMOEP'SON  (1). 

c 

C  CmLLLO  OYJ  call  HT-ITOl  (!4,M,  X,PHI,3f  A,P> 

c 

C  N  IS  THZ  TCIMI  «  ?  AMpL  LS 

C  M  IS  THE  Of  D-R  PLUS  Ot  b 

C  X()  *••;£  TMr  t MPU7  VALUES 

C  PHIO  ARE  r  M  ;  CALC  UL  'TEC  t  U7  CCOKRCL*  Ti  Q»*  OFFriciENTS 
C  BT A (  )  ACE  TH  =  CALC!JLAT£0  t  ST  I  FAT  ION  E^OP  70"rFIClENTS 
C  PI)  APS  THE  T  MNCVA  T I C  M  VAr.IA  NCI  <~'OWIO  V*L'JFC 

C 

C  1.  Ai-.GERSO'i,  •).  ."CALFULATIOf  CF  FIL1 Eh  C  OE  FFTC  x  ENTS  FOR 
C  'IAXI-J-t  -•MTROP’J  S  =  CCTRAL  nf  ALYSIS",  Sf  OPHYS.,  VOL. 39, 

C  M  L:  •  1 ,  > 3  *■ '  2 , F  E  0  7  A  • 

C 

4  4*  4~  *  *  *  •  *  *♦#44*  -  -  •*  !•**■<'  *  J-  ♦•i  4  **  »  *•*.■*•'  *■♦***•' '*4^44***»»  ** 

c 

31  4ENSI0N  X  (1)  ,PHI(l)  ,?T  A  (1)  fP(l)  ,31  (Tir  )  ,B2(r0O>  ,A?7A(  li?) 
C  FIND  P(l)  Afn  PHI  ( 1 ) 

SUM= 

00  IP  1=1,  N 

10  SIJM=$'JM  +  X  (I ) *  *  2 • 

P<1)  =S'JM/ FLOAT  (N) 

PHI(1)=P(1) 

3T  A  ( 1 )  =  1  • 

IFCt.ET.l)  RETURN 
C  NOW  complete  the  REGRESSION 

MM  1=  M- 1 
DO  2c  J  =  1 , NM1 
**1  (1)  =X  (I  > 

23  82  Ci ) =X (I +1 ) 

1=2 

30  SM0M  =  .  , 

soen=:. 

Nr1IPl=S-I  +  1 
00  L  *i  J  =  1,NMIP1 
SHOM=SMOM  *3 1  ( J)  *  n  2  (  J) 

40  S3EM=33ENMni(  J)  ■  '  2  .♦  P2(  J)*-  2.) 

5TA(I)=2. ‘SNOrt/SCEN 
P(I)=P(T*1)“  (1,-PTA  (I  )■*•■  2.) 

SUH=C. 

oo  -;g  <=?, i 

50  iUM=S'JM-PHI  «I-<>  1  >  ♦  3TA(K) 

PHI(I)=«?U1 
iF(il,LT,  2)  RETURN 
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j»r  U.E0.2>G0  TO  ■*  •„ 

T*11=I-t 

00  ...  J  K=?,IMl 

6u  9TA(^)  =  i9TA(K)-37A  ( I)  '  A^T  A  (I  -K*l> 

if  (i  .Eo.M)  return 
7(1  1*1*1 

iMtrl-l 

oo  oT 

83  A* TA«)=',-A  (<) 

NMIPi*N-I+l 

00  i">  *=1,NM1P1 

01(K)  =  ni  m  -ATT A  (I-D-  Sr  (K) 

90  32(<)=°?(<  +  l)-i3*  A(T-l)  «1(<*1) 

GO  TO  ? 

end 


Subroutine  PAPO 

Given  a  sample  of  an  input  time  process,  this  subroutine 

determines  an  estimate  of  the  frequency  spectrum  using  the 

techniques  developed  by  Papoulis  Eqs(56)  through  (60). 

.*  *  •  •  »  ***'•»■  .  1. **♦•*♦*»»* 

c 

SUBROUTINE  FAP0(^,17tL.Fr  ,  FHI  ,DT,  X  *  P°,  FRM  ,  ;  JT ) 

C 

C  PAPO  CALC"l:tts  the  SPECTnAL  CEuSlTY  USING  THE  METHODS 
C  0E3C rvi  BED  EY  °AP0ULIS  (i). 

C 

C  CALL  £0  BY  i  PALL  PAPG  (N  ,1  T  ,C  Eiv,  FHI , 01  ,  X,  !» , FRO,  ITT) 

C 

CM  IS  THE  NUMBER  THAT  IS  USED  TO  CALCULATE  "HE  TOTAL 

C  N'J~»BEP  OF  SAMPLES.  N  =  c"“  M 

C  IT  IS  THE  TIME  SAMPLE  HINrCW  *0R  THE  ITTERA'ION 
C  DER  IS  RE  01 1 1'1  ED  MINIMUM  CHANGE  14  «5E  ONE  ITTEkATI^N 

C  TO  THE  NEXT 

C  FHI  IS  THI  HT  GHC  ST  F^OUcNCY  OF  xNT  £R  ”S” 

C  DT  IS  THE  TIME  SAMPLE  INTERVAL 

C  XO  ARE  THE  r  I  ME  SAMPLE  iNrl'i  VARIABLES 
C  PPO  ARE  T H£  CALCULATED  POWER  SPECTkUM  OUTPUT  VALUES 
C  FROO  ARE  THF  FREQUENCY  SAMPLE  POTMTS 

C  ITT  IS  THE  TOTAL  NUKBt*  OF  ITTSkATiONS  *0  MEET  THt  OER 
C  VALUE,  OR  the  PIFfJLT  Value  it j 

C 

C  REQUIRED  ROUTINES*  FFT 

c 

C  1.  PAPOULIS,  A .  ,'*A  NEW  ALGOMTHM  IN  SPECTRAL  ANALYSIS  .MD  3AN0- 
.  C  LiMlTEO  '‘XTRAP0LATILNm,1F  Et  TRANS.  CIP.  SYS .,  VOL  .  CAS-  22, 

C  NC.9,S£P  fr . 

C 

c 
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Oi SESSION  FPO(i>  ,*<i)  ,FF  (1) 

complex  p’jn (?;:), ft 

N=2  ’  M 
nO  2P1=  l+N/2 

C  GENERATE  FLUENCY  FOIST 3 
CF*(2.'FHi) /FLOAT  (N ) 

F;>n=. . 

0<>  b>  I*l,w 
FPO(I)=FREO 
5  FREO=FPFn+OF 

C  GEi4Z-aTE  aSITTal  SAMPLE  WINCGW, 
C  FWS( i) = ' 

DG  i«i  T  =  1,N 
IF(a.TT.I~)GC  TO  2. 

Fwsm=*  m 

Gu  TO  1  : 

2G  Fu*j(;)*r.i  t.> 
in  COST  T*l  JF 

C  GET  THE  FFT 

CALL  crT(FWU,Mf i.) 

Ii  T  =  1 

3G  DC  J*NS2P1,N 
4  0  FWN  ( J)  =  (  C  .  »  t« ) 

CALL  FFT (FWN,M,-1.) 

C  TEST  FOR  M5E  LESS  THAW  PEL 


+■ 

* 


& 
v  v 


<*y 


s 


FSN(I)=X(I), (* .LE.IT),  ELSE 


cK  *  -•  • 

TK1=«j. 

00  5  L  =  1 ,1T 

A  =  CONJG(FMS  (L) ) 

Ai*A *  P  MS ( L ) 

A2s<AP3  <Y(L>>  v2.) 
aV/AL=A1*«2-2.-  30.TT(iii)-?PF.T{  A2) 
EFsER+CMl+AVAL)  /2. 

45  1  S  l=  ?  V  A  L 

GLTER=E?M1-A3i (EF ) 

Er-sABS  (PP) 

ERN1=ER 

IFdTT.LE.lJGO  TO  -7 
C  IF  MSE  LESS  THAN  OER,  DOSE 
lFOLTER.LE.OEs)  GC  TO  L? 

C  IF  1-.:  ITTE“  ATI  ONS,  OGME 
Ic  (ITT.  3c.l  ■;«-■?)  GO  TO 
C  ELSE,  DO  ANOTHER  ITEFATION 
47  00  31  J=i  » IT 

50  FWS(J)sX(J) 

CALL  FFT (PWK,H,i. ) 
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ITT*ITT  ♦1 

3 0  TO  ?•» 

C  OOrtE*  GET  ^o^OUENCY  E.P EC T  »dJ  P 
60  i,',L'-  rrT  (CWN»  M»  1  •  ) 

!)£•**£* 

00  i  i  r=i,N 

pp(I)srw»HI)*CONJG(FWK(J  )  ) 
70  DD(*)=iO-’ (PP(I) ) 

KETU^O 

EN3 


Subroutine  PDF 

PDF  generate  uniform  distribution  deviates  using  the 
following  algorithm: 

INPUT  Y 
OUTPUT  F 
A-0. 

B-l. 

1.  If  Y  is  greater  than  A,  GO  to  4 

2.  F*0 

3.  GO  TO  8 
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*  c*r* 


c 

SUBROUTINE  Pne(Y,F) 

C 

C  PDF  COMPUTE?  the  THEOlvETI  CAL  UNIFORM  DIE  ~RI  **U”  ION  FUNCTION 
C 

C  CALLcD  PY?  FILL  PDF(Y,F) 

c 

C  Y  IS  INPUT 
C  F  IS  OUTPUT 

c 

c 

J  • 

5=1. 

ifiy.gt.dgc  to  •■ 

F=  "i . 

GO  TO  1  > 

5  iF (Y .Li*  5) GO  TO  1 . 

r=l. 

GO  TO  1" 

13  F=(Y-.1)/(P-m) 

15  RETURN 

ENO 

Subroutine  PSPEC 

This  subroutine  calculates  the  estimated  power  spectrum 
using  Eq( 18)  from  the  values  of  §n  and  PM+1  determined  in 
the  routine  MEMCOl. 


P(o>) 


PM+lAt 


M 


ne 

n*l 


-ju>nAt 


(18) 


where  At  is  the  input  sampling  interval  and  M  is  the  AR 
model  order. 
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. . 


far* 


c 

sut-.o'j t  r  »•  s  psp zc  ( : f i s ,  :t  :  r » fl  ) ,  ^'T ,  3,i , "t  a , P3 ,  f  > 

c 

C  PS3EC  COM=flTFS  THE  ES7IHA1SL'  FQW t \  S’ESTPU*  iriNG  THE 
C  MEN  EQUATION  (1). 

C 

C  CALL  EC  9Y»  CALL  PSP  7 C ( I P 1 S, J CKO t JT ♦ FL C, FHT , c V , 9Tfi , PS , F> 
C 

C  IPTS  IS  ?H z  M'MflFR  3^  PflHlJ  COMPUTE! 

C  IO*J  IS  TUI  Hf.-Ocl  O'  £  Ck  PLUi  ONE 
C  OT  IS  THE  T  T  H  ^  SAMPLE  ll\T  tRVAL 
C  FLO  iS  TH"  LOW  FREQUENCY 
C  FHI  13  TH-  Miocsi  FREPUENCY 

C  PH  IS  THE  ET^Ck  1NNOVMIOUS  V^uVi CL  (POWER)  pOR 
C  FuR  OR"’?  m 

C  8TA()  A«E  THT  £$1 I HAT i JN  CCr  FFICIENI  S 
C  PSD  a  RE  THE  CjmPUtEP  °QW  if.  SPEC  T\U-  POT*T? 

C  FC)  ARE  THE  F'ECUEHCY  PClH'lS 

C 

C  1.  A  i»OE;S  jN  i  N.  ,**C  AL  CULAT  xGt  Lc  FILTER  CCEcPTr  1  ENTS  FO’’ 
C  MfiXI'HV,  EM7R0PV  SprCTRPL  ANAl YaI3"»  jf  PPHY? . »  VOL • 39  » 

C  .MG.  l,  :i-’2,FEF  7v  . 

c 
c 

0IMEN3I  ON  8TA(i)  ,PS  (1)  ,F  (1) 

COMPLEX  A,1 
OF=(FHr-rLO)/IPTS 
PI  =  3.1  •»l';«r2tL33«98 
FnEO=F.O 
00  1C  t=i,ifts 
F(  I)  =rr  Er 

IF(I0PO.G7.1)G0  TO  ! 

PatDsPH'  OT 
GO  TO  1"' 

5  X=2. 4  PI  *  F ( 1 )  *  OT 

csiJN= : . 
ssum=:. 

00  2G  J=2,I0RD 

CS0Ms:3'J-+9Ti  (  J)  COS  (  (J-1)*X) 

20  SMlJM>5fifv+8TA  (  J)  SI  M  (  (J-l  DX) 
csuH=i.-rsuM 
AsCMPLxrr'-UH.ssu1  ) 

0=CONJG ( A ) 

A0VAL=a*  ° 

PS(I)=(PH  0T)/A3VAL 
15  PRE0=F5Er»*0F 

10  CONTIUJF 

RETURN 
ENO 


/ 


/ 

<*■  (  Cj 

&  ? 


,•»  ;• 

Q  y> 

•S  a 

*4  ^ 

5* 
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Appendix  C 


Additional  Plots  from  Analysis 
of  Known  Process  in  Section  II 

The  following  figures  are  additional  analysis  results 
from  examples  two  and  three  of  Section  II.  Figures  41 
through  47  show  additional  results  from  analysis  of  the  ARMA 
process.  The  additional  data  from  the  sum  of  two  cosines  is 
provided  in  Figures  48  through  54. 


estimation  of  system 

ORDER  US  I  NO  flKIUKE 
FINAL  PREDICTOR  ERR 
(FfE) 
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Figure  41.  Amplitude  of  FPE  as  a  Function  of  Model  Order  for  ARMA  Model 
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MEM  Technique,  ARMA  Model 
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Figure  45.  Estimated  Mean  of  Estimated  Spectrum  Using  the 
Papoulis  Technique,  ARMA  Model 
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Figure  48.  Amplitude  of  FPE  as  a  Function  of  Model  Order  for 

Sum  of  Sinusoids  Input 
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Figure  50.  Probability  of  Kolmogorov-Smirnov  Two  Sided  Test  as 
Function  of  Model  Order  for  Sum  of  Sinusoids  Input 
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Figure  51.  Estimated  Mean  of  Estimated  Spectrum  Using  the 
MEM  Technique,  Sum  of  Sinusoids  Input 


Figure  53.  Estimated  Variance  of  Estimated  Spectrum  Using  the 
MEM  Technique,  Sum  of  Sinusoids  Input 
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Figure  54.  Estimated  Variance  of  Estimated  Spectrum  Using  the 
Papoulis  Technique,  Sum  of  Sinusoids  Input 
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